include "/clima/model64/togradsnew.lib/tograds.inc" character * 8 Cdate8 character * 3 cend,cihr dimension s(8),pradi8(8),pdens8(8),id7(7),ar7(7),qext(8) real, ALLOCATABLE, dimension(:,:,:,:):: s8, varpl real, ALLOCATABLE, dimension(:,:,:) :: varsfc, z, dl, sfc, dd, dw real, ALLOCATABLE, dimension(:,:,:) :: pint, sold, zold real, ALLOCATABLE, dimension(:,:,:) :: qoqi, ice, cwm, dns real, ALLOCATABLE, dimension(:,:,:) :: tznd, tsrf real, ALLOCATABLE, dimension(:,:,:) :: u, v, q, t,omg real, ALLOCATABLE, dimension(:,:,:) :: fimm,fimma,fimmb,ext real, ALLOCATABLE, dimension(:,:) :: z5,z7,icel,dmsk,psfc real, ALLOCATABLE, dimension(:,:) :: aot,sfcsvi,ddsvi,dwsvi,hgt real, ALLOCATABLE, dimension(:,:) :: t2m,psl,lsm real, ALLOCATABLE, dimension(:,:) :: sfc4, dlsvi, dl4, u10, v10 real, ALLOCATABLE, dimension(:,:) :: acp, fimml, fimmal,fimmbl real,parameter :: und=-999., und1=9.999e+20 DATA PRADI8/0.15,0.25,0.47,0.80,1.36,2.29,3.93,7.24/ DATA PDENS8/2.50,2.50,2.50,2.50,2.65,2.65,2.65,2.65/ DATA QEXT /1.536,2.816,3.086,2.583,2.345,2.240,2.174,2.110/ integer:: nlevmx,it double precision, ALLOCATABLE, dimension(:,:,:) :: fix double precision:: ff real:: raot,s8e(8) data id7/09,05,25,00,0,0,0/ real:: dl8(8) wctln(1:1) = ' ' call getarg(1,cdate8) call getarg(2,cend) call getarg(3,cihr) read (cend,*) iend read (cihr ,*) inchr read (Cdate8(1:2),'(i2)') id7(1) read (Cdate8(3:4),'(i2)') id7(2) ! id7(3) = 1 read (Cdate8(5:6),'(i2)') id7(3) read (Cdate8(7:8),'(i2)') id7(4) print*, "radim posshort za ", id7(1:3) c print *,' id7=',id7 call datetohr(id7(3),id7(2),id7(1),id7(4),id7(5),IHRS00) call datetohr(id7(3),id7(2),id7(1),id7(4),id7(5),IHRS) call RGRADS (id7,81, 001, 0, 1, 1, 1, 0., DDD, NLRET) IMU = rnx_ctl !(1) JMU = rny_ctl !(1) NZ = nz_ctl (1) ALLOCATE(u (IMU,JMU,NZ)) ALLOCATE(v (IMU,JMU,NZ)) ALLOCATE(z (IMU,JMU,NZ)) ALLOCATE(t (IMU,JMU,NZ)) ALLOCATE(q (IMU,JMU,NZ)) ALLOCATE(omg (IMU,JMU,NZ)) ALLOCATE(cwm (IMU,JMU,NZ)) ALLOCATE(ice (IMU,JMU,NZ)) ALLOCATE(qoqi (IMU,JMU,NZ)) ALLOCATE(pint (IMU,JMU,NZ)) ALLOCATE(tznd (IMU,JMU,NZ)) ALLOCATE(tsrf (IMU,JMU,NZ)) ALLOCATE(fimm (IMU,JMU,NZ)) ALLOCATE(fimma (IMU,JMU,NZ)) ALLOCATE(fimmb (IMU,JMU,NZ)) ALLOCATE(dns (IMU,JMU,NZ)) ALLOCATE(ext (IMU,JMU,NZ)) ALLOCATE(fix (IMU,JMU,NZ)) ALLOCATE(zold (IMU,JMU,NZ+1)) ALLOCATE(sold (IMU,JMU,NZ+1)) ALLOCATE(s8 (IMU,JMU,NZ,8)) ALLOCATE(dl (IMU,JMU,8)) ALLOCATE(sfc (IMU,JMU,8)) ALLOCATE(dd (IMU,JMU,8)) ALLOCATE(dw (IMU,JMU,8)) ALLOCATE(icel (IMU,JMU)) ALLOCATE(psfc (IMU,JMU)) ALLOCATE(aot (IMU,JMU)) ALLOCATE(hgt (IMU,JMU)) ALLOCATE(sfcsvi(IMU,JMU)) ALLOCATE(ddsvi (IMU,JMU)) ALLOCATE(dwsvi (IMU,JMU)) ALLOCATE(dlsvi (IMU,JMU)) ALLOCATE(sfc4 (IMU,JMU)) ALLOCATE(dl4 (IMU,JMU)) ALLOCATE(u10 (IMU,JMU)) ALLOCATE(v10 (IMU,JMU)) ALLOCATE(acp (IMU,JMU)) ALLOCATE(z5 (IMU,JMU)) ALLOCATE(z7 (IMU,JMU)) ALLOCATE(t2m (IMU,JMU)) ALLOCATE(dmsk (IMU,JMU)) ALLOCATE(fimml (IMU,JMU)) ALLOCATE(fimmal(IMU,JMU)) ALLOCATE(fimmbl(IMU,JMU)) ALLOCATE(psl (IMU,JMU)) ALLOCATE(lsm (IMU,JMU)) ! NOB_CTL(NCTLMX),XOB_CTL(NOBSMX),YOB_CTL(NOBSMX) ALLOCATE(varpl (IMU,JMU,NZ,16)) ALLOCATE(varsfc(IMU,JMU,16)) ! fimml=und ! icel=und ! dmsk=und ! dlsvi=und dlmdU = rxi_ctl !(1) dphdu = ryi_ctl !(1) bwU = rx0_ctl ! (1) bsU = ry0_ctl ! (1) beU = bwU + (IMU-1)*dlmdU bnU = bsU + (JMU-1)*dphdU print *,imu,"f",bwu,beu,dlmdu print *,jmu,"f",bsu,bnu,dphdu print *,NOB_CTL(1) do k=1,NOB_CTL(1) print *,XOB_CTL(k),YOB_CTL(k) enddo ! ------ Reykjavik --------------------------------------- jbg = (64.13-bsU)/dphdu+1 ibg = (-21.87-bwU)/dlmdu+1 rlat = bsU+jbg*dphdu rlon = bwU+ibg*dlmdu ! print*, 'Reykjavik RLON RLAT IBG JBG :', rlon, rlat,ibg,jbg ! open(23,file="Reykjavik_sfcconc.txt") ! -------------------------------------------------------- open (30,file="posshort.gdat",status='unknown' + ,form="unformatted",access="direct",recl=IMU*JMU) ! open (60,file="posshort_sfc.gdat",status='unknown' ! + ,form="unformatted",access="direct",recl=IMU*JMU) ! ! open (70,file="posshort_vis.gdat",status='unknown' ! + ,form="unformatted",access="direct",recl=IMU*JMU) c bojan do 500 IHT=IHRS+3,IHRS+IEND,inchr !****************** GLAVNA PETLJA PO VREMENU ******************************** do 500 IHT=IHRS,IHRS+IEND,inchr !**************************************************************************** JHT = IHT call hrtodate(JHT,id7(5),0,id7(3),id7(2),id7(1),id7(4)) write(334,*) 'VALID: ', id7(1:5) call RGRADS (id7, 81, 1, 0,IMU,JMU,1, 0., lsm, NLRET) call RGRADS (id7, 07, 1, 0,IMU,JMU,1, 0., hgt, NLRET) call RGRADS (id7, 07,100,500,IMU,JMU,1,500., z5, NLRET) call RGRADS (id7, 07,100,700,IMU,JMU,1,700., z7, NLRET) call RGRADS (id7,200, 1, 0,IMU,JMU,1, 0., dmsk, NLRET) call RGRADS (id7, 1, 1, 0,IMU,JMU,1, 0., psfc, NLRET) call RGRADS (id7, 2,102, 0,IMU,JMU,1, 0., psl, NLRET) do l=1, NZ call RGRADS (id7, 39,109, l,IMU,JMU,1,1.*l, omg(:,:,l), NLRET) call RGRADS (id7, 7,109, l,IMU,JMU,1,1.*l, z(:,:,l), NLRET) call RGRADS (id7, 11,109, l,IMU,JMU,1,1.*l, t(:,:,l), NLRET) call RGRADS (id7, 51,109, l,IMU,JMU,1,1.*l, q(:,:,l), NLRET) call RGRADS (id7,153,109, l,IMU,JMU,1,1.*l, cwm(:,:,l), NLRET) call RGRADS (id7, 1,109, l,IMU,JMU,1,1.*l, pint(:,:,l), NLRET) ! call RGRADS (id7, 58,109, l,IMU,JMU,1,1.*l, ice(:,:,l), NLRET) call RGRADS (id7, 33,109, l,IMU,JMU,1,1.*l, u(:,:,l), NLRET) call RGRADS (id7, 34,109, l,IMU,JMU,1,1.*l, v(:,:,l), NLRET) do k=1,8 call RGRADS (id7,200+k,109,l,IMU,JMU,1,1.*l, s8(:,:,l,k),NLRET) enddo ! k where(s8(:,:,:,:).gt.0.9*und1)s8=0. enddo ! lev cdl idl=221-1 isf=211-1 idd=231-1 idw=241-1 do k=1,8 call RGRADS(id7,idl+k, 1, 0,IMU,JMU,1, 0., dl(:,:,k), NLRET) call RGRADS(id7,isf+k, 1, 0,IMU,JMU,1, 0., sfc(:,:,k), NLRET) call RGRADS(id7,idd+k, 1, 0,IMU,JMU,1, 0., dd(:,:,k), NLRET) call RGRADS(id7,idw+k, 1, 0,IMU,JMU,1, 0., dw(:,:,k), NLRET) enddo call RGRADS(id7, 33,105,10,IMU,JMU,1,10., u10, NLRET) call RGRADS(id7, 34,105,10,IMU,JMU,1,10., v10, NLRET) call RGRADS(id7, 61, 1, 0,IMU,JMU,1, 0., acp, NLRET) call RGRADS(id7, 11,105, 2,IMU,JMU,1, 2., t2m, NLRET) ! do j=1,jmu ! do i=1,imu ! rxx=dl(i,j,1) ! if(rxx.ne.0.and.rxx.lt.9.99E+19) then ! write(333,*) i,j,dl(i,j,1) ! endif ! enddo ! enddo ! stop "CCCCCCCCCCCCCCCCCCCCCCCCC" !====== otkoci ako zelis proveru za neko polje =========== ! print*,'provera za t' ! do j=1,jmu ! do i=1,imu ! do l=1,nz ! rxx=t(i,j,l) ! if(q(i,j,l).ne.0) then ! ! if(rxx.ne.0.and.rxx.lt.9.99E+19) then ! ! !print*, i,j,l,q(i,j,l) ! write(333,*) i,j,l,rxx ! endif ! enddo ! enddo ! enddo !========================================================= do j=1,jmu do i=1,imu do l=1,nz zold(i,j,l) = -hgt(i,j)+z(i,j,l) sold(i,j,l) = 0 do k=1, 8 sold(i,j,l) = sold(i,j,l) +s8(i,j,l,k) enddo enddo enddo enddo do j=1,jmu do i=1,imu dl(i,j,:) = 0. icel(i,j) = 0. zb = hgt(i,j) do l=1,nz-1 c dz = z(i,j,l)-z0 zt=0.5*(z(i,j,l)+z(i,j,l+1)) dz = zt-zb pp=0.5*(PINT(i,j,l)+PINT(i,j,l+1)) tt=t(i,j,l) rho=pp/(tt*287.) !bojan PAZI OVDE VEROVATNO TREBA DA SE MNOZI SA GUSTINOM icel(i,j)=icel(i,j)+dz*ice(i,j,l)/rho ! ovo je masa leda po m2 (icel) do k=1,8 dl(i,j,k)=dl(i,j,k)+dz*s8(i,j,l,k) enddo ! k zb = zt enddo enddo ! i enddo ! j aot = 0. sfcsvi = 0. ddsvi = 0. dwsvi = 0. do j=1,jmu do i=1,imu if( hgt(i,j) .lt.5000.and.hgt(i,j).ge.-200 ) then c call checkp(imu,jmu,u10,-100.,100.,und) ! call dl2aotxx(8, dl(i,j,1:8),aot(i,j)) dl8(1:8)=dl(i,j,1:8) ! call dl2aotxx(8, dl(i,j,1:8),raot) call dl2aotxx(8, dl8(1:8),raot) aot(i,j)=raot sfcsvi(i,j) = sum(sfc(i,j,1:8)) dlsvi (i,j) = sum(dl (i,j,1:8)) ddsvi (i,j) = sum(dd (i,j,1:8)) dwsvi (i,j) = sum(dw (i,j,1:8)) sfc4 (i,j) = sum(sfc(i,j,1:4)) dl4 (i,j) = sum(dl (i,j,1:4)) endif ! if(sfcsvi(i,j)*1e9.lt.0.0000001) sfcsvi(i,j)=0. ! if(i==ibg.and.j==jbg) then ! write(23,445)id7(1:4),sfcsvi(i,j)*1e9 ! print*, sfcsvi(i,j)*1e9 ! stop ! endif end do end do 445 format(I4,2x,3(I2.2,2x),2x,f6.1) do j=1,jmu do i=1,imu ! --------------------------------- FIMML (i,j) = 0. FIMMaL(i,j) = 0. FIMMbL(i,j) = 0. zb = hgt(i,j) pb = psfc(i,j) do L=1,NZ-1 zt=0.5*(z(i,j,l)+z(i,j,l+1)) dz = zt-zb ! --------------------------------- TSRF (i,j,l) = 0. TZND (i,j,l) = 0. ! TZND(:,:,:) = 0. FIMM (i,j,l) = 0. FIMMa(i,j,l) = 0. FIMMb(i,j,l) = 0. DNS (i,j,l) = 0. qoqi (i,j,l) = 0. ! --------------------------------- DO N=1, 8 !7 ! -------------------------------------------------------------- S(N)=S8(I,J,L,N) ! *1E9 !converted to ug/m^3 !!! write(334,*)i,j,l,n,s8(i,j,l,n) ! N particles per unite volume ! ZND(I,J,L,N)=(S(N)*1.E-6) ! */(PDENS8(N)*1.E+6*4./3.*3.14*(PRADI8(N)*1.E-6)**3) !#/m^3 ! -------------------------------------------------------------- ZND=(S(N)*1.) * /(PDENS8(N)*1.E+3*4./3.*3.14*(PRADI8(N)*1.E-6)**3) !#/m^3 TZND(I,J,L)=TZND(I,J,L)+ZND !#/m^3 ! -------------------------------------------------------------- ! surface concentration for bins SRF=ZND*4*3.14*(PRADI8(N)*1E-6)**2 !m^2/m^3 ! total N of particles per unite volume ! -------------------------------------------------------------- ! total surface concentration TSRF(I,J,L)=TSRF(I,J,L)+SRF !m^2/m^3 ! -------------------------------------------------------------- ENDDO ! N ! write(333,*)i,j,l,tznd(i,j,l)*1e-6 TT=T(I,J,L)-273.16 ! pp=PINT(i,j,l) pp=0.5*(PINT(i,j,l)+pb) ! if(q(i,j,l).lt.9.99E+19) print*,'qq',i,j,l, q(i,j,l) CALL tp2qw(pp,tt,qw,qi,qint) qoqi(i,j,l) = q(i,j,l)/qi ! if(it==1) then ! write(776,*) i,j,l,q(i,j,l),qi ! endif ! !#################################################################################### IF(q(i,j,l).ge.0.8*qi)THEN ! IF(q(i,j,l).ge.1.1*qi)THEN ! IF(TT.LT.-5..AND.TT.GT.-35.)THEN IF(TT.LT.-15..AND.TT.GT.-35.)THEN ! ===================== NIMAND ================================ DNS(I,J,L)=EXP(-0.517*TT+8.934) ! [1/m^2] FIMM(I,J,L)=DNS(I,J,L)*TSRF(I,J,L)/1000 ! #/litres ! ============================================================== ! ===================== De Mott A ============================ ! alternative - modified DeMott (convert tznd to #/cm^3)!FIMM2 in #/litres FIMMa(I,J,L)=0.0008*(10**(-0.2*(TT+9.7))) & *(TZND(I,J,L)/1000000.)**1.25 ! ============================================================== ! ===================== De Mott B ============================ bb=1.25 gg=0.46 ddd=-11.6 FIMMb(I,J,L)=3.*((TZND(I,J,L)/1000000.)**bb)*exp(gg*(-TT)+ddd) ! ============================================================== ! if(l==21) then ! if(FIMM2(I,J,L).gt.0.)then ! r1=FIMM2(I,J,21) ! endif ! endif ! if(l==21)write(555,*)i,j,r1-FIMM2(I,J,l) ! -------------------------------------------------------------- ELSE IF(TT.LT.-35..AND.TT.GT.-55.)THEN ! -------------------------------------------------------------- ! Steinke 2014 si=1.03 si=qoqi(i,j,l) xth=-TT+(si-1)*100 ! xth=-1.085*TT+0.815*(si-1)*100 ! alternativno DNS(I,J,L)=1.88*1.e5*EXP(0.2659*xth) FIMM(I,J,L)=DNS(I,J,L)*TSRF(I,J,L)/1000. ! #/litres FIMMa(I,J,L)=FIMM(I,J,L) !!!! Demot+St OBJEDINJENO DEMOTT I STEINKE !!! FIMMb(I,J,L)=FIMM(I,J,L) !!!! Demot+St OBJEDINJENO DEMOTT I STEINKE !!! ! FIMM2(I,J,L)=0.0008*(10**(-0.2*(TT+9.7))) ! demot ! &*(TZND(I,J,L)/1000000.)**1.25 ENDIF ENDIF FIMML(i,j) = FIMML(i,j) + 1000.* dz * FIMM(I,J,L) ! #/m2 LOAD FIMMaL(i,j) = FIMMaL(i,j) + 1000.*dz * FIMMa(I,J,L) ! #/m2 LOAD FIMMbL(i,j) = FIMMbL(i,j) + 1000.*dz * FIMMb(I,J,L) ! #/m2 LOAD ! write(999,*)fimm(i,j,l) ! if(fimmb(i,j,l).gt.1e16) then ! write(992,*)fimmb(i,j,l), TZND(I,J,L), DNS(i,j,L), TT ! endif ! if(omg(i,j,l).lt.0.)then it=(IHT-IHRS)/inchr+1 !(iyyv-yyst)*12+immv zb = zt pb = pint(i,j,l) s8e(1:8)=s8(i,j,l,1:8) ! call extinc(8, s8(i,j,l,1:8),ext(i,j,l)) call extinc(8, s8e(1:8),ext(i,j,l)) fix(i,j,l)=-(omg(i,j,l)**3)*log10(fimm(i,j,l)+1) if(omg(i,j,l)>=0) fix(i,j,l)=0. enddo ! L ! if(fimmbl(i,j).gt.1e16) then ! ! write(993,553)i,j,log10(fimmbl(i,j)+1), TZND(I,J,15:21),TT ! write(993,*)i,j,TZND(I,J,15:15) ! endif ! !553 format(I3.3,3x,I3.3,3x,f15.8,2x,1f11.8,4x,f7.3) ! 553 format(I3.3,3x,I3.3,3x,f7.3,7(f12.7,2x),f7.3) enddo ! i enddo ! j ! stop "BBBBBBBBBBBBBBBB" !================ provera za fimm =================== ! open (77,file="fimm.gdat" ! + ,form="unformatted" ! + ,access="direct",recl=IMU*JMU) ! ! do L=NL1,NL2 ! it1=(IHT-IHRS00)/inchr ! irec1=it1+L-8+1 ! write(77,rec=irec1) fimm(:,:,l) ! ! write(666,*)it1,l,irec1 !,minval(fimm),maxval(fimm) ! enddo !====================================================== !====================================================== ! upis u pos.gdat !------------------------------------------------------ ! irec=(iens-1)*nt*nvar*nlev+(it-1)*nlev*nvar+(iv-1)*nlev+il ! it=(iyyv-yyst)*12+immv ! irec=(it-1)*nvarisfc+(it-1)*nvarpl*nlev+iv ! und1 je nedefinisana vrednost u pos.ctl where(abs(sold).gt.1e15) sold=und1 ! where(fix.lt.0.) fix=0. c open (30,file="posshortsvifimm.gdat" ! open (30,file="posshort.gdat" ! + ,form="unformatted",access="direct",recl=IMU*JMU) it=(IHT-IHRS00)/inchr+1 print*, it nvar =22 nvarsfc =11 nvarpl =11 nlev =28 varsfc(:,:,1) =lsm varsfc(:,:,2) =sfcsvi varsfc(:,:,3) =dlsvi varsfc(:,:,4) =fimmbl varsfc(:,:,5) =t2m varsfc(:,:,6) =acp varsfc(:,:,7) =aot varsfc(:,:,8) =u10 varsfc(:,:,9) =v10 varsfc(:,:,10) =ddsvi varsfc(:,:,11) =dwsvi varpl(:,:,:,1) =u varpl(:,:,:,2) =v varpl(:,:,:,3) =q varpl(:,:,:,4) =t varpl(:,:,:,5) =omg varpl(:,:,:,6) =z varpl(:,:,:,7) =zold varpl(:,:,:,8) =sold varpl(:,:,:,9) =fimmb varpl(:,:,:,10)=tznd varpl(:,:,:,11)=fix do iv=1,nvar do il=1,nlev if(iv.le.nvarsfc) then !======================= SURFACE ============================= irec=(it-1)*nvarsfc+(it-1)*nvarpl*nlev+iv irecw=irec write (30,rec=irecw) varsfc(:,:,iv) !============================================================= else !======================== LEVELS ============================= irec=(it-1)*nvarsfc+(it-1)*nvarpl*nlev+ & nvarsfc+(iv-nvarsfc-1)*nlev+il irecw=irec write(30,rec=irecw) varpl(:,:,il,iv-nvarsfc) !============================================================= endif enddo !lev enddo !var !!upis samo prizemna polja ! do iv=1,nvarsfc ! irecs=(it-1)*nvarsfc+iv ! write(60,rec=irecs)varsfc(:,:,iv) ! enddo !!upis samo visinska polja ! do iv=1,nvarpl ! do il=1,nlev ! irecl=(it-1)*nlev*nvarpl+(iv-1)*nlev+il ! write(70,rec=irecl) varpl(:,:,il,iv) ! enddo ! enddo ! 500 continue end !111 format(i2,1x,i2,1x,i2,1x,i2,3x,i2.2,2x,i2,1x,i2,1x,i4.4) C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine checkp(im,jm,v,bl,bu,und) dimension v(im,jm) do j=1,jm do i=1,im if( v(i,j).ge.bl.and.v(i,j).lt.bu) then else v(i,j) = und endif enddo enddo return end c================================================================== subroutine dl2aotxx(KPS,dload,raot) dimension QEXT550(8), PRADI(8),PDENS(8),dload(KPS) DATA QEXT550 /1.373,3.303,3.245,2.413,2.262,2.260,2.162,2.108/ !cbojan ovo je pre bilo DATA PRADI /0.15,0.25,0.45,0.78,1.32,2.24,3.80,7.11/ !microns effective DATA PRADI /0.15,0.25,0.47,0.80,1.36,2.29,3.93,7.24/ DATA PDENS /2.50,2.50,2.50,2.50,2.65,2.65,2.65,2.65/ !g/cm3 aot8 = 0 do n =1, 8 spart = dload(n) aot8=aot8 + + spart * 1000*3*QEXT550(N)/(4*PRADI(N)*PDENS(N)) enddo raot = aot8 c apradi apdens aqext 2.012500 2.575000 2.390750 c print *,k,1000*spart,OTHIC550 return end subroutine extinc(KPS,stot,ext) dimension QEXT(8), PRADI(8),PDENS(8),stot(KPS) DATA QEXT /1.536,2.816,3.086,2.583,2.345,2.240,2.174,2.110/ !Tegenlacis DATA PRADI /0.15,0.25,0.47,0.80,1.36,2.29,3.93,7.24/ DATA PDENS /2.50,2.50,2.50,2.50,2.65,2.65,2.65,2.65/ !g/cm3 ext = 0 do n =1, 8 spart = stot(n) ext=ext + + spart * 1000*3*QEXT(N)/(4*PRADI(N)*PDENS(N)) enddo ! aot = aot8 c apradi apdens aqext 2.012500 2.575000 2.390750 c print *,k,1000*spart,OTHIC550 return end !call spline(lold,zslh,vij,y2,lm,zuv,vtil,pp,qq) subroutine spline(nold,xold,yold,y2,nnew,xnew,ynew,p,q) ! ! ****************************************************************** ! * * ! * this is a one-dimensional cubic spline fitting routine * ! * programed for a small scalar machine. * ! * * ! * programer: z. janjic, yugoslav fed. Hydromet. Inst., beograd * ! * * ! * * ! * * ! * nold - number of given values of the function. Must be ge 3. * ! * xold - locations of the points at which the values of the * ! * function are given. Must be in ascending order. * ! * yold - the given values of the function at the points xold. * ! * y2 - the second derivatives at the points xold. If natural * ! * spline is fitted y2(1)=0. And y2(nold)=0. Must be * ! * specified. * ! * nnew - number of values of the function to be calculated. * ! * xnew - locations of the points at which the values of the * ! * function are calculated. Xnew(k) must be ge xold(1) * ! * and le xold(nold). * ! * ynew - the values of the function to be calculated. * ! * p, q - auxiliary vectors of the length nold-2. * ! * * ! ****************************************************************** ! d i m e n s i o n 2 xold(nold),yold(nold),y2(nold),p(nold),q(nold) 3,xnew(nnew),ynew(nnew) !----------------------------------------------------------------------- noldm1=nold-1 ! dxl=xold(2)-xold(1) dxr=xold(3)-xold(2) dydxl=(yold(2)-yold(1))/dxl dydxr=(yold(3)-yold(2))/dxr rtdxc=.5/(dxl+dxr) ! p(1)= rtdxc*(6.*(dydxr-dydxl)-dxl*y2(1)) q(1)=-rtdxc*dxr ! if(nold.eq.3) go to 700 !----------------------------------------------------------------------- k=3 ! 100 dxl=dxr dydxl=dydxr dxr=xold(k+1)-xold(k) dydxr=(yold(k+1)-yold(k))/dxr dxc=dxl+dxr den=1./(dxl*q(k-2)+dxc+dxc) ! p(k-1)= den*(6.*(dydxr-dydxl)-dxl*p(k-2)) q(k-1)=-den*dxr ! k=k+1 if(k.lt.nold) go to 100 !----------------------------------------------------------------------- 700 k=noldm1 ! 200 y2(k)=p(k-1)+q(k-1)*y2(k+1) ! k=k-1 if(k.gt.1) go to 200 !----------------------------------------------------------------------- k1=1 ! 300 xk=xnew(k1) ! do 400 k2=2,nold if(xold(k2).le.xk) go to 400 kold=k2-1 go to 450 400 continue ynew(k1)=yold(nold) go to 600 450 if(k1.eq.1) go to 500 if(k.eq.kold) go to 550 ! 500 k=kold ! y2k=y2(k) y2kp1=y2(k+1) dx=xold(k+1)-xold(k) rdx=1./dx ! ak=.1666667*rdx*(y2kp1-y2k) bk=.5*y2k ck=rdx*(yold(k+1)-yold(k))-.1666667*dx*(y2kp1+y2k+y2k) ! 550 x=xk-xold(k) xsq=x*x ! ynew(k1)=ak*xsq*x+bk*xsq+ck*x+yold(k) ! 600 k1=k1+1 if(k1.le.nnew) go to 300 !----------------------------------------------------------------------- return end subroutine tp2qw(pp,tcel,qw,qi,qint) C*********************************************************************** P A R A M E T E R & (H1=1.E0,H2=2.E0 &, D00=0.E0,D125=.125E0,D5=0.5E0 &, A1=610.78E0,A2=17.2693882E0,A3=273.16E0,A4=35.86E0 &, PQ0=379.90516E0,TRESH=.95E0 &, CP=1004.6E0,ELWV=2.50E6,ELIV=2.834E6,ROW=1.E3,G=9.8E0 &, EPSQ=2.E-12,DLDT=2274.E0,TM10=263.16E0,R=287.0E0 &, CPR=CP*R,RCPR=H1/(CPR)) P A R A M E T E R & (ARCP=A2*(A3-A4)/CP,RCP=H1/CP,PQ0C=PQ0*TRESH,RROG=H1/(ROW*G)) C-------------------- t=273.16+tcel TMT0=T-273.16 TMT15=AMIN1(TMT0,-15.) AI=0.008855 BI=1. IF(TMT0.LT.-20.)THEN AI=0.007225 BI=0.9674 ENDIF QW=PQ0/PP*EXP(A2*(T-A3)/(T-A4)) QI=QW*(BI+AI*AMIN1(TMT0,0.)) QINT=QW*(1.-0.00032*TMT15*(TMT15+15.)) IF(TMT0.LE.-40.)QINT=QI return end