calculate_watercycle.f90 7.24 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 Oct 13, 2017 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 37 `````` !***************************************************************************** ! * ! 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 integer :: itime,i,j,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 Nov 29, 2017 39 `````` integer :: interp_time `````` Sabine committed Oct 13, 2017 40 41 `````` integer :: partnumber,forparticle ! -1 inititalize q, partnumber real do deposition real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,topo `````` Sabine committed Jan 17, 2018 42 `````` real centerposx_real, centerposy_real, centerposx, centerposy, diff, xl(2), yl(2), prec `````` Sabine committed Oct 12, 2017 43 44 45 46 47 48 49 50 51 52 `````` real qv1(2),qvprof(2),qvi,dz,dz1,dz2 ! Some variables needed for temporal interpolation !************************************************* dt1=real(itime-memtime(1)) dt2=real(memtime(2)-itime) dtt=1./(dt1+dt2) `````` Sabine committed Nov 29, 2017 53 54 55 56 57 `````` indexh=memind(2) if (abs(memtime(1)-itime).lt.abs(memtime(2)-itime)) & indexh=memind(1) `````` Sabine committed Jan 17, 2018 58 `````` ! if called from timemanager a loop over all particle is necessary `````` Sabine committed Oct 16, 2017 59 60 61 62 `````` ! if particle is released, the initialization has just to be done ! for the single particle !****************************************************************** `````` Sabine committed Oct 13, 2017 63 64 `````` if (partnumber.eq.-1) then forparticle=numpart `````` Sabine committed Oct 16, 2017 65 `````` write (*,*) 'calculate watercycle, ',numpart,partnumber,forparticle,itime `````` Sabine committed Oct 13, 2017 66 67 68 `````` else forparticle=1 endif `````` Sabine committed Oct 12, 2017 69 `````` `````` Sabine committed Oct 16, 2017 70 `````` ! Loop about all particles `````` Sabine committed Oct 12, 2017 71 `````` !************************* `````` Sabine committed Oct 13, 2017 72 73 74 75 76 77 `````` do j=1,forparticle if (partnumber.eq.-1) then i=j ! loop over all particle else i=partnumber ! just for one particle `````` Sabine committed Feb 07, 2018 78 `````` firsttimestep(i)=.true. `````` Sabine committed Oct 13, 2017 79 `````` endif `````` Sabine committed Oct 12, 2017 80 `````` `````` Sabine committed Oct 16, 2017 81 82 `````` ! Take only valid particles, itime invalid when particle released! !***************************************************************** `````` Sabine committed Oct 12, 2017 83 `````` `````` Sabine committed Nov 29, 2017 84 `````` if ((itra1(i).eq.itime).or.(forparticle.eq.1)) then `````` Sabine committed Oct 12, 2017 85 86 87 88 89 `````` !***************************************************************************** ! Interpolate several variables (PV, specific humidity, etc.) to particle position !***************************************************************************** `````` Sabine committed Oct 16, 2017 90 91 92 93 94 95 96 97 ``````! Determine age class of the particle - nage is used for the kernel !****************************************************************** itage=abs(itra1(i)-itramem(i)) do nage=1,nageclass if (itage.lt.lage(nage)) goto 33 end do 33 continue `````` Sabine committed Oct 12, 2017 98 99 100 101 `````` ix=xtra1(i) jy=ytra1(i) ixp=ix+1 jyp=jy+1 `````` Sabine committed Oct 16, 2017 102 `````` if (jy.eq.numygrid) jyp=numygrid `````` Sabine committed Oct 12, 2017 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 `````` ddx=xtra1(i)-real(ix) ddy=ytra1(i)-real(jy) rddx=1.-ddx rddy=1.-ddy p1=rddx*rddy p2=ddx*rddy p3=rddx*ddy p4=ddx*ddy ! Topography !*********** topo=p1*oro(ix ,jy) & + p2*oro(ixp,jy) & + p3*oro(ix ,jyp) & + p4*oro(ixp,jyp) ! Potential vorticity, specific humidity, temperature, and density !***************************************************************** 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) ! Specific humidity qv1(m)=p1*qv(ix ,jy ,ind,indexh) & +p2*qv(ixp,jy ,ind,indexh) & +p3*qv(ix ,jyp,ind,indexh) & +p4*qv(ixp,jyp,ind,indexh) 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 Nov 29, 2017 153 154 155 156 157 158 159 `````` prec=p1*(lsprec(ix ,jy ,1,indexh)+convprec(ix ,jy ,1,indexh)) & +p2*(lsprec(ix ,jy ,1,indexh)+convprec(ix ,jy ,1,indexh)) & +p3*(lsprec(ix ,jy ,1,indexh)+convprec(ix ,jy ,1,indexh)) & +p4*(lsprec(ix ,jy ,1,indexh)+convprec(ix ,jy ,1,indexh)) `````` Sabine committed Jan 24, 2018 160 `````` if (partnumber.eq.-1) then ! not the initialisation stage `````` Sabine committed Oct 12, 2017 161 162 163 164 `````` ! Calculate humidity differences !******************************* `````` Sabine committed Jan 24, 2018 165 `````` diff=(qvi-val_q(i))*xmass1(i,1)*ldirect `````` Sabine committed Oct 12, 2017 166 167 168 169 `````` ! Determine center of mass position on earth and average height !************************************************************** `````` Sabine committed Jan 24, 2018 170 171 172 173 174 `````` 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 endif `````` Sabine committed Oct 12, 2017 175 `````` `````` Sabine committed Oct 13, 2017 176 `````` ! save the old values - for particle after release - particlenumber = -1 `````` Sabine committed Oct 12, 2017 177 178 `````` !******************** `````` Sabine committed Jan 24, 2018 179 180 181 `````` xtra1_q(i)=xtra1(i) ytra1_q(i)=ytra1(i) val_q(i)=qvi `````` Sabine committed Oct 12, 2017 182 `````` `````` Sabine committed Nov 29, 2017 183 184 `````` if (partnumber.eq.-1) then ! not the initialisation stage `````` Sabine committed Feb 02, 2018 185 186 ``````! Do this only for the 2nd timestep i ! there has to be a deltaq which is negative and the column has to have precipitation `````` Sabine committed Feb 07, 2018 187 188 189 ``````! e_minus_p units is mm/day .. convert to mm/h if ( ((diff.ge.0).or.(e_minus_p(ix,jy)/24.gt.-2.0)) .and. firsttimestep(i)) then ! if ( ((diff.ge.0)) .and. firsttimestep(i)) then `````` Sabine committed Jan 24, 2018 190 `````` itra1(i)=-999999999 `````` Sabine committed Feb 07, 2018 191 192 193 194 195 196 197 198 199 200 201 202 203 ``````! write(*,*) 'Trajectory terminated: ',i,diff,e_minus_p(ix,jy),firsttimestep(i) else ! write(*,*) 'Trajectory accounted: ',i,diff,e_minus_p(ix,jy),firsttimestep(i),& ! val_q(i),qvi,xmass1(i,1) prec_q(i)=prec call centerofmass(xl,yl,2,centerposx_real,centerposy_real) centerposx=(centerposx_real-xlon0)/dx centerposy=(centerposy_real-ylat0)/dy if (diff.gt.0) then call drydepokernel(1,diff,centerposx,centerposy,nage,1) else call wetdepokernel(1,abs(diff),centerposx,centerposy,nage,1) endif `````` Sabine committed Jan 24, 2018 204 `````` endif `````` Sabine committed Feb 07, 2018 205 `````` firsttimestep(i)=.false. `````` Sabine committed Nov 29, 2017 206 `````` endif `````` Sabine committed Oct 16, 2017 207 `````` `````` Sabine committed Nov 29, 2017 208 `````` endif `````` Sabine committed Oct 12, 2017 209 210 211 212 `````` end do end subroutine calculate_watercycle``````