c hmm.f by Y. ashie NILIM-MLIT (2010.12.5) c •Û…«ŒšÞ‚Ì”M…•ª“¯ŽžˆÚ“®‰ð̓vƒƒOƒ‰ƒ€ c 1•ª–ˆ‚ɉ·“xAŠÜ…—¦A…•ªƒ|ƒeƒ“ƒVƒƒƒ‹‚ðo—ÍB c@@ ‰ŠúðŒF‰·“xŠÜ…—¦‹ÏˆêAã–ÊF’f”M’fŽ¼A’ê–ÊF’f”M‹z… c @-----------------•Ï”à–¾-------------------------- c tp(i)F‘æiƒƒbƒVƒ…‚̉·“x[K] c@@ wp(i)F‘æiƒƒbƒVƒ…‚Ì‘ÌÏŠÜ…—¦[-] c tp1(i)FVŽž‚Ì‘æiƒƒbƒVƒ…‚̉·“x[K] c@@ wp1(i)FVŽž‚Ì‘æiƒƒbƒVƒ…‚Ì‘ÌÏŠÜ…—¦[-] real*8 tp(100),tp1(100),wp(100),wp1(100),pffsave(100) open(2,file='hmm_1.txt') ppi=3.1415926535 eps=0.95 sumconv=0. sumevap=0. sumcond=0. sumbeta=0. ii=0 xle=600.*1000.*4.18 xlsample=30./1000. nsample=31 dx=xlsample/(float(nsample-2)+.5) deltt=.05 wmax=0.568 xx=60./deltt jj=int(xx) write(*,*)'jj deltt',jj,deltt tme=0. tmemax=1.*24.*3600.*5. do 10 i=1,nsample tp(i)=30. wp(i)=0.01 10 continue c--- sabun ------ 40 ii=ii+1 tme=deltt*float(ii) if (tme.ge.tmemax+deltt*1.2) goto 1000 if (mod(ii,jj).eq.1) then ta=30. xa=.02 xi=0. xl=0. vel=1. endif alfc=7.2*vel+4.2 alfx=alfc/1000. do 20 i=2,nsample-1 ww=wp(i-1) tt=tp(i-1) call calw(ww,tt,pff,rhh,dww,dtt,dwvv,dtvv,xkk,rmm,crr) pff1=pff dw1=dww dt1=dtt dwv1=dwvv dtv1=dtvv xk1=xkk rm1=rmm cr1=crr ww=wp(i) tt=tp(i) call calw(ww,tt,pff,rhh,dww,dtt,dwvv,dtvv,xkk,rmm,crr) pff0=pff dw0=dww dt0=dtt dwv0=dwvv dtv0=dtvv xk0=xkk rm0=rmm cr0=crr ww=wp(i+1) tt=tp(i+1) call calw(ww,tt,pff,rhh,dww,dtt,dwvv,dtvv,xkk,rmm,crr) pff2=pff dw2=dww dt2=dtt dwv2=dwvv dtv2=dtvv xk2=xkk rm2=rmm cr2=crr wp1(i)=wp(i)+deltt/1000./dx* 1 ((dw2+dw0+dwv2+dwv0)/2.*(wp(i+1)-wp(i))/dx 1 +(dw1+dw0+dwv1+dwv0)/2.*(wp(i-1)-wp(i))/dx 1 +(dt2+dt0+dtv2+dtv0)/2.*(tp(i+1)-tp(i))/dx 1 +(dt1+dt0+dtv1+dtv0)/2.*(tp(i-1)-tp(i))/dx 1 +9.80665*xk1-9.80665*xk2) tp1(i)=tp(i)+((rm2+rm0)/2.*(tp(i+1)-tp(i)) 1 +(dwv2+dwv0)/2.*(wp(i+1)-wp(i))*xle 1 +(dwv1+dwv0)/2.*(wp(i-1)-wp(i))*xle 1 +(rm1+rm0)/2.*(tp(i-1)-tp(i)))*deltt/cr0/dx/dx if (i.eq.nsample-1) then wp1(i)=wp(i)+deltt/1000./dx* 1 (0.*(dw2+dw0+dwv2+dwv0)/2.*(wp(i+1)-wp(i))/dx 1 +(dw1+dw0+dwv1+dwv0)/2.*(wp(i-1)-wp(i))/dx 1 +0.*(dt2+dt0+dtv2+dtv0)/2.*(tp(i+1)-tp(i))/dx 1 +(dt1+dt0+dtv1+dtv0)/2.*(tp(i-1)-tp(i))/dx 1 +9.80665*xk1-9.80665*xk2*0.) tp1(i)=tp(i)+(0.*(rm2+rm0)/2.*(tp(i+1)-tp(i)) 1 +0.*(dwv2+dwv0)/2.*(wp(i+1)-wp(i))*xle 1 +(dwv1+dwv0)/2.*(wp(i-1)-wp(i))*xle 1 +(rm1+rm0)/2.*(tp(i-1)-tp(i)))*deltt/cr0/dx/dx endif if (wp1(i).gt.0.5675) wp1(i)=.5675 if (wp1(i).lt..0021) wp1(i)=.0021 if (tp1(i).gt.100.) tp1(i)=100. if (tp1(i).lt.0.) tp1(i)=0. 20 continue c i=1 ww=wp(1) tt=tp(1) call calw(ww,tt,pff,rhh,dww,dtt,dwvv,dtvv,xkk,rmm,crr) pff0=pff dw0=dww dt0=dtt dwv0=dwvv dtv0=dtvv xk0=xkk rm0=rmm cr0=crr xss1=6.11*10.**(7.5*tp(1)/(tp(1)+237.3)) xss=0.622*xss1/(1013.25-xss1) xss2=6.11*10.**(7.5*tp(1)/(tp(1)+237.3))*rhh xs=0.622*xss2/(1013.25-xss2) beta=(xs-xa)/(xss-xa) if (beta.lt.0.) beta=0. if (beta.gt.1.) beta=1. if (wp(1).lt.0.005) beta=0. ww=wp(2) tt=tp(2) call calw(ww,tt,pff,rhh,dww,dtt,dwvv,dtvv,xkk,rmm,crr) pff2=pff dw2=dww dt2=dtt dwv2=dwvv dtv2=dtvv xk2=xkk rm2=rmm cr2=crr r=0. d=0. cond=(rm2+rm0)/2.*(tp(2)-tp(1))/dx 1 +(dtv2+dtv0)/2.*(wp(2)-wp(1))/dx*xle evap=alfx*(xa-xss)*beta conv=alfc*(ta-tp(1)) evap=0. cond=0. conv=0. sumconv=sumconv+conv sumevap=sumevap+evap sumcond=sumcond+cond sumbeta=sumbeta+beta wp1(1)=wp(1)+(r-d+evap 1 +(dw2+dw0+dwv2+dwv0)/2.*(wp(2)-wp(1))/dx 1 +(dt2+dt0+dtv2+dtv0)/2.*(tp(2)-tp(1))/dx 1 -9.8055*xk2)/dx*deltt/(1000./2.) tp1(1)=tp(1)+(0.8*xi+eps*xl-eps*4.88e-8/.86*(tp(1)+273.15)**4 1 +cond+evap*xle+conv)*deltt/dx/(cr0/2.) if (wp1(1).gt..5675) wp1(1)=.5675 if (wp1(1).lt.0.002) wp1(1)=.0021 if (tp1(1).gt.100.) tp1(1)=100. if (tp1(1).lt.0.) tp1(1)=0. c-----------------(4/4)-------------------------- if (tme.gt.0.) then do 25 i=nsample-1,nsample wp1(i)=.5 tp1(i)=30. 25 continue endif c----------------------------------------------- do 30 i=1,nsample tp(i)=tp1(i) wp(i)=wp1(i) ww=wp(i) tt=tp(i) call calw(ww,tt,pff,rhh,dww,dtt,dwvv,dtvv,xkk,rmm,crr) pffsave(i)=pff 30 continue if (mod(ii,jj).eq.0) then if (tme.ge.0.) then do 50 i=1,nsample 50 continue if (mod(ii,24*60*60).eq.0) 1 write(*,*)'I=',ii,' Day=',tme/60./60./24.,deltt,tme sumconv=sumconv/float(jj) sumevap=sumevap/float(jj) sumcond=sumcond/float(jj) sumbeta=sumbeta/float(jj) write(2,100)tme,ta,tp1(1),tp1(2),tp1(3),tp1(4),tp1(5), 1 tp1(10),tp1(15),tp1(20),tp1(25),tp1(30),tp1(31), 1 wp1(1),wp1(2),wp1(3),wp1(4),wp1(5), 1 wp1(6),wp1(7),wp1(8),wp1(9),wp1(10),wp1(11), 1 wp1(12),wp1(13),wp1(14),wp1(15),wp1(16),wp1(17), 1 wp1(18),wp1(19),wp1(20),wp1(21),wp1(22),wp1(23), 1 wp1(24),wp1(25),wp1(26),wp1(27),wp1(28),wp1(29), 1 wp1(30),wp1(31),pffsave(1),pffsave(2), 1 pffsave(3),pffsave(4),pffsave(5), 1 pffsave(6),pffsave(7),pffsave(8), 1 pffsave(9),pffsave(10),pffsave(11), 1 pffsave(12),pffsave(13),pffsave(14), 1 pffsave(15),pffsave(16),pffsave(17), 1 pffsave(18),pffsave(19),pffsave(20), 1 pffsave(21),pffsave(22),pffsave(23), 1 pffsave(24),pffsave(25),pffsave(26), 1 pffsave(27),pffsave(28),pffsave(29), 1 pffsave(30),pffsave(31) 1 ,0.8*xi+eps*xl-0.*eps*4.88e-8/.86*(tp(1)+273.15)**4 1 ,-sumconv,-sumevap*xle,sumcond,sumbeta sumconv=0. sumevap=0. sumcond=0. sumbeta=0. 100 format(200e14.5) endif endif goto 40 1000 stop end subroutine calw(ww,tt,pff,rhh,dww,dtt,dwvv,dtvv,xkk,rmm,crr) c @-----------------•Ï”à–¾-------------------------- c wwF‘ÌÏŠÜ…—¦[-] c ttF‰·“x[K] c@@ pffF…•ªƒ|ƒeƒ“ƒVƒƒƒ‹[J/kg] c rhh:‘Š‘ÎŽ¼“x[-] c@@ dww:ŠÜ…—¦Œù”z‚̉t…“`“±Œn”[kg/ms] c@@ dtt:‰·“xŒù”z‚̉t…“`“±Œn”[kg/msK] c@@ dwvv:ŠÜ…—¦Œù”z‚Ì…ö‹C“`“±Œn”[kg/ms] c@@ dtvv:‰·“xŒù”z‚Ì…ö‹C“`“±Œn”[kg/msK] c@ @xkk:“§…ŒW”[kgs/m3] c rmm:”M“`“±—¦[W/mK] c crr:‘ÌÏ”M—e—Ê[J/m3K] if (ww.gt.0.5675) ww=0.5675 if (ww.lt..0021) ww=0.0021 if (tt.gt.100.) tt=100. if (tt.lt.0.) tt=0. ss=(ww-0.002)/(0.568-0.002) rmm=0.006*ww*100.+0.14 crr=41800.*ww*100.+548250. pff= (ss**(-1./0.23)-1.)**(1./1.3)/11.83 xkk= ss**0.5*(1.-(1.-ss**(1./0.23))**0.23)**2*3.1/9.80665 dww=xkk*(1./0.23/1.3/11.83/(0.568-0.002)*(ss**(-1./0.23)-1.) 1 **(1./1.3-1)*ss**(-1./0.23-1.)) dtt=0. dwvv=0. dtvv=0. pff=-pff rhh=exp(pff/(461.52)/(273.15+tt)) return end