calculate_watercycle.f90 10.2 KB
 Sabine committed Oct 12, 2017 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 ``````!********************************************************************** ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 * ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * ! * ! This file is part of FLEXPART. * ! * ! FLEXPART is free software: you can redistribute it and/or modify * ! it under the terms of the GNU General Public License as published by* ! the Free Software Foundation, either version 3 of the License, or * ! (at your option) any later version. * ! * ! FLEXPART is distributed in the hope that it will be useful, * ! but WITHOUT ANY WARRANTY; without even the implied warranty of * ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * ! GNU General Public License for more details. * ! * ! You should have received a copy of the GNU General Public License * ! along with FLEXPART. If not, see . * !********************************************************************** `````` Sabine committed Aug 23, 2018 22 ``````subroutine calculate_watercycle(partnumber,itime) `````` Sabine committed Oct 12, 2017 23 24 25 26 27 28 29 30 31 32 33 34 35 36 `````` !***************************************************************************** ! * ! calculate difference in q and dump it to the dry/wetdepofield * ! VAR: real, allocatable, dimension(:) :: xtra1_q, ytra1_q, val_q * ! * !***************************************************************************** use par_mod use com_mod implicit none real(kind=dp) :: jul `````` Sabine committed Jun 06, 2018 37 `````` integer :: itime,i,j,k,jjjjmmdd,ihmmss,numshortout,numshortall `````` Sabine committed Oct 16, 2017 38 `````` integer :: ix,jy,ixp,jyp,ind,indz,indzp,il,m,indexh,itage,nage `````` Sabine committed Jun 06, 2018 39 `````` integer :: interp_time, ngrid `````` Sabine committed Jun 30, 2018 40 `````` integer :: partnumber,forparticle,totalparticle ! -1 inititalize q, partnumber real do deposition `````` Sabine committed Jun 06, 2018 41 `````` real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,xtn,ytn `````` Sabine committed Feb 08, 2018 42 `````` real centerposx_real, centerposy_real, centerposx, centerposy, diff, xl(2), yl(2) `````` Sabine committed Oct 12, 2017 43 `````` real qv1(2),qvprof(2),qvi,dz,dz1,dz2 `````` Sabine committed Oct 30, 2018 44 `````` real e_minus_p1(2), e_minus_pi, threshold `````` Sabine committed Jun 30, 2018 45 `````` `````` Sabine committed Oct 30, 2018 46 `````` real precip_reference(360,180),eminusp_saved `````` Sabine committed Jul 10, 2018 47 `````` integer ixr,jyr,idumx,idumy,istep, watersynctime `````` Sabine committed Jun 30, 2018 48 `````` `````` Sabine committed Oct 12, 2017 49 `````` `````` Sabine committed Oct 30, 2018 50 51 52 `````` watersynctime=3600*3 ! was loutstep before, but for monthly run not possible ! watersynctime=3600 ! watersynctime=loutstep `````` Sabine committed Apr 19, 2019 53 `````` threshold = -.5 ! units mm/h `````` Sabine committed Sep 26, 2018 54 `````` `````` Sabine committed Apr 19, 2019 55 `````` if ((partnumber.eq.1).and.(itime.eq.0)) write(*,*) 'Watersynctime ',watersynctime,lsynctime, ' s, threshold: ',threshold, 'mm/h' `````` Sabine committed Aug 23, 2018 56 `````` `````` Sabine committed Oct 12, 2017 57 58 59 60 61 62 63 `````` ! Some variables needed for temporal interpolation !************************************************* dt1=real(itime-memtime(1)) dt2=real(memtime(2)-itime) dtt=1./(dt1+dt2) `````` Sabine committed Jan 17, 2018 64 `````` ! if called from timemanager a loop over all particle is necessary `````` Sabine committed Oct 16, 2017 65 66 67 68 `````` ! if particle is released, the initialization has just to be done ! for the single particle !****************************************************************** `````` Sabine committed Oct 13, 2017 69 70 `````` if (partnumber.eq.-1) then forparticle=numpart `````` Sabine committed Jun 30, 2018 71 `````` totalparticle=0 `````` Sabine committed Jul 10, 2018 72 `````` istep=int(itime/watersynctime) `````` Sabine committed Oct 30, 2018 73 ``````! write(*,*) 'Model step:',istep `````` Sabine committed Jun 30, 2018 74 75 76 77 78 79 80 81 82 83 84 85 86 87 ``````! open(43,file='precip_ref_input.dat',status='old') ! write(*,*) 'Reading eminusp from file' ! do 75 jyr=1,90 ! do 75 ixr=1,360 !75 read(43,'(2i4,e15.3)') idumx,idumy,precip_reference(ixr,jyr) ! close(43) ! open(43,file='precip_ref_plus_eminusp.dat') ! write(*,*) 'Writing eminusp from file' ! do 76 jyr=1,90 ! do 76 ixr=1,360 !6 write(43,'(2i4,3e15.3)') ixr,jyr,precip_reference(ixr,jyr), & ! e_minus_p(ixr,jyr,1), e_minus_p(ixr,jyr+90,1) ! close(43) ! write(*,*) 'Written eminusp from file' `````` Sabine committed Oct 13, 2017 88 89 90 `````` else forparticle=1 endif `````` Sabine committed Oct 12, 2017 91 `````` `````` Sabine committed Feb 08, 2018 92 `````` ! Loop about all particles `````` Sabine committed Oct 12, 2017 93 `````` !************************* `````` Sabine committed Oct 13, 2017 94 95 96 97 98 99 100 `````` do j=1,forparticle if (partnumber.eq.-1) then i=j ! loop over all particle else i=partnumber ! just for one particle endif `````` Sabine committed Oct 12, 2017 101 `````` `````` Sabine committed Jul 06, 2018 102 ``````! write(*,*) 'particle: ',i,xtra1_q(i),xtra1(i),itime,itra1(i),itramem(i),status_q(i) `````` Sabine committed Oct 16, 2017 103 104 `````` ! Take only valid particles, itime invalid when particle released! !***************************************************************** `````` Sabine committed Feb 08, 2018 105 `````` if (itra1(i).eq.itime) then ! trajectory itra1 is not -999999 `````` Sabine committed Oct 12, 2017 106 `````` `````` Sabine committed Feb 08, 2018 107 108 `````` ! Determine age class of the particle - nage is used for the kernel !****************************************************************** `````` Sabine committed Oct 16, 2017 109 110 111 112 113 `````` itage=abs(itra1(i)-itramem(i)) do nage=1,nageclass if (itage.lt.lage(nage)) goto 33 end do 33 continue `````` Sabine committed Jul 06, 2018 114 `````` `````` Sabine committed Jun 30, 2018 115 `````` `````` Sabine committed Jun 06, 2018 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 `````` ! Determine which nesting level to be used !***************************************** ngrid=0 do k=numbnests,1,-1 if ((xtra1(i).gt.xln(k)).and.(xtra1(i).lt.xrn(k)).and. & (ytra1(i).gt.yln(k)).and.(ytra1(i).lt.yrn(k))) then ngrid=k goto 23 endif end do 23 continue ! Determine nested grid coordinates !********************************** if (ngrid.gt.0) then xtn=(xtra1(i)-xln(ngrid))*xresoln(ngrid) ytn=(ytra1(i)-yln(ngrid))*yresoln(ngrid) ix=int(xtn) jy=int(ytn) ddy=ytn-real(jy) ddx=xtn-real(ix) else `````` Sabine committed Jun 11, 2018 140 `````` ix=int(xtra1(i)) `````` Sabine committed Jun 06, 2018 141 142 143 144 145 `````` jy=int(ytra1(i)) ddx=xtra1(i)-real(ix) ddy=ytra1(i)-real(jy) endif `````` Sabine committed Oct 12, 2017 146 147 `````` ixp=ix+1 jyp=jy+1 `````` Sabine committed Oct 16, 2017 148 `````` if (jy.eq.numygrid) jyp=numygrid `````` Sabine committed Oct 12, 2017 149 150 151 152 153 154 155 `````` rddx=1.-ddx rddy=1.-ddy p1=rddx*rddy p2=ddx*rddy p3=rddx*ddy p4=ddx*ddy `````` Sabine committed Jun 30, 2018 156 157 `````` ! Interpolate specific humidity - code from partoutput.f90 !********************************************************* `````` Sabine committed Oct 12, 2017 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 `````` do il=2,nz if (height(il).gt.ztra1(i)) then indz=il-1 indzp=il goto 6 endif end do 6 continue dz1=ztra1(i)-height(indz) dz2=height(indzp)-ztra1(i) dz=1./(dz1+dz2) do ind=indz,indzp do m=1,2 indexh=memind(m) `````` Sabine committed Jun 06, 2018 177 178 `````` if (ngrid.eq.0) then qv1(m)=p1*qv(ix ,jy ,ind,indexh) & `````` Sabine committed Oct 12, 2017 179 180 181 `````` +p2*qv(ixp,jy ,ind,indexh) & +p3*qv(ix ,jyp,ind,indexh) & +p4*qv(ixp,jyp,ind,indexh) `````` Sabine committed Jun 06, 2018 182 183 184 185 186 187 `````` else qv1(m)=p1*qvn(ix ,jy ,ind,indexh,ngrid) & +p2*qvn(ixp,jy ,ind,indexh,ngrid) & +p3*qvn(ix ,jyp,ind,indexh,ngrid) & +p4*qvn(ixp,jyp,ind,indexh,ngrid) endif `````` Sabine committed Oct 12, 2017 188 189 190 191 192 193 `````` end do qvprof(ind-indz+1)=(qv1(1)*dt2+qv1(2)*dt1)*dtt end do qvi=(dz1*qvprof(2)+dz2*qvprof(1))*dz `````` Sabine committed Jul 10, 2018 194 `````` if ( (mod(itage,watersynctime).eq.0).and.(itage.gt.0) ) then !trajectory is loutstep after simulation start, not a new particle `````` Sabine committed Jul 06, 2018 195 `````` `````` Sabine committed Feb 12, 2018 196 197 198 `````` do m=1,2 indexh=memind(m) `````` Sabine committed Jun 06, 2018 199 200 `````` if (ngrid.eq.0) then e_minus_p1(m)=p1*e_minus_p(ix ,jy ,indexh) & `````` Sabine committed Feb 12, 2018 201 202 203 `````` +p2*e_minus_p(ixp,jy ,indexh) & +p3*e_minus_p(ix ,jyp,indexh) & +p4*e_minus_p(ixp,jyp,indexh) `````` Sabine committed Jun 06, 2018 204 205 206 207 208 209 `````` else e_minus_p1(m)=p1*e_minus_p_nest(ix ,jy ,indexh,ngrid) & +p2*e_minus_p_nest(ixp,jy ,indexh,ngrid) & +p3*e_minus_p_nest(ix ,jyp,indexh,ngrid) & +p4*e_minus_p_nest(ixp,jyp,indexh,ngrid) endif `````` Sabine committed Feb 12, 2018 210 211 `````` end do e_minus_pi=(e_minus_p1(1)*dt2+dt1*e_minus_p1(2))*dtt `````` Sabine committed Jun 30, 2018 212 `````` eminusp_saved=e_minus_pi `````` Sabine committed Feb 12, 2018 213 `````` `````` Sabine committed Feb 08, 2018 214 215 `````` ! Calculate humidity differences !******************************* `````` Sabine committed Nov 29, 2017 216 `````` `````` Sabine committed Jul 06, 2018 217 `````` diff=(qvi-val_q(i))*xmass1(i,1)*ldirect `````` Sabine committed Oct 12, 2017 218 219 220 221 `````` ! Determine center of mass position on earth and average height !************************************************************** `````` Sabine committed Feb 08, 2018 222 223 224 225 `````` xl(1)=xlon0+xtra1_q(i)*dx xl(2)=xlon0+xtra1(i)*dx yl(1)=ylat0+ytra1_q(i)*dy yl(2)=ylat0+ytra1(i)*dy `````` Sabine committed Oct 12, 2017 226 `````` `````` Sabine committed Feb 08, 2018 227 ``````! Do this only for the 2nd timestep `````` Sabine committed Feb 02, 2018 228 ``````! there has to be a deltaq which is negative and the column has to have precipitation `````` Sabine committed Feb 07, 2018 229 ``````! e_minus_p units is mm/day .. convert to mm/h `````` Sabine committed Feb 08, 2018 230 `````` `````` Sabine committed Jun 30, 2018 231 232 233 ``````! e_minus_pi=precip_reference(ix+1,jy-90+2)*(-24.) ! First timestep and criteria fullfilled, keep it or throw it? `````` Sabine committed Oct 30, 2018 234 `````` if ( ((diff.ge.0).or.((e_minus_pi/24).gt.threshold)) .and. (status_q(i).eq.1) ) then `````` Sabine committed Aug 23, 2018 235 ``````! if ( ((diff.ge.0).or.((e_minus_pi/24).gt.-2.0)) .and. (status_q(i).eq.1) ) then `````` Sabine committed Aug 21, 2018 236 ``````! if ( ((diff.ge.0)) .and. (status_q(i).eq.1) ) then `````` Sabine committed Jul 06, 2018 237 238 239 240 241 ``````! if (status_q(i).eq.-100) then ! never! all trajectories. ! write(*,*) 'terminated: ',i,diff,e_minus_pi,itage,itramem(i),status_q(i),qvi,val_q(i),xmass1(i,1) status_q(i)=-9 ! this particle is invalid itra1(i)=-999999999 else `````` Sabine committed Jun 30, 2018 242 243 244 245 246 247 `````` totalparticle=totalparticle+1 call centerofmass(xl,yl,2,centerposx_real,centerposy_real) centerposx=(centerposx_real-xlon0)/dx centerposy=(centerposy_real-ylat0)/dy `````` Sabine committed Jul 06, 2018 248 ``````! write(*,*) 'accounted: ',i,itime,itra1(i),itage,diff,itramem(i),status_q(i) `````` Sabine committed Jun 30, 2018 249 250 251 252 `````` if (diff.gt.0) then call drydepokernel(1,diff,centerposx,centerposy,nage,1) if (nested_output.eq.1) & call drydepokernel_nest(1,diff,centerposx,centerposy,nage,1) `````` Sabine committed Jul 06, 2018 253 `````` else `````` Sabine committed Jun 30, 2018 254 255 256 `````` call wetdepokernel(1,abs(diff),centerposx,centerposy,nage,1) if (nested_output.eq.1) & call wetdepokernel_nest(1,abs(diff),centerposx,centerposy,nage,1) `````` Sabine committed Jul 06, 2018 257 `````` endif `````` Sabine committed Jul 10, 2018 258 259 260 261 ``````! if (istep.lt.100) then ! waterfieldp(istep,ix+1,jy+1)=waterfieldp(istep,ix+1,jy+1)+xmass1(i,1) ! endif ! waterfielde(ix+1,jy+1)=waterfielde(ix+1,jy+1)+xmass1(i,1) `````` Sabine committed Jul 06, 2018 262 263 264 265 266 `````` xtra1_q(i)=xtra1(i) ytra1_q(i)=ytra1(i) val_q(i)=qvi status_q(i)=2 endif `````` Sabine committed Jun 30, 2018 267 268 `````` endif ! at louttimestep interval `````` Sabine committed Jul 06, 2018 269 270 271 272 273 274 275 276 277 `````` if (status_q(i).eq.-1) then ! first time, a value is saved ! write(*,*) 'saved: ',i,xtra1(i),ytra1(i),qvi,itime,itage,itramem(i),status_q(i) xtra1_q(i)=xtra1(i) ytra1_q(i)=ytra1(i) val_q(i)=qvi status_q(i)=1 endif `````` Sabine committed Jun 30, 2018 278 `````` endif ! for valid trajectories - itra1.eq.itime `````` Sabine committed Jul 06, 2018 279 `````` `````` Sabine committed Oct 12, 2017 280 `````` end do `````` Sabine committed Jul 06, 2018 281 `````` if ((partnumber.eq.-1).and.(totalparticle.gt.0)) then `````` Sabine committed Oct 30, 2018 282 ``````! write (*,*) 'calculated watercycle, ',totalparticle,numpart,partnumber,forparticle,itime `````` Sabine committed Jun 30, 2018 283 `````` endif `````` Sabine committed Oct 12, 2017 284 285 286 `````` end subroutine calculate_watercycle``````