added also euler_nest

src/euler_rain_nests.f90

+      subroutine euler_rain_nest(itime)
+!*******************************************************************************
+! Compute E-P between 2 time steps from moisture budget in an atmospheric column
+! 2 Terms: 1. Moisture column tendency between the 2 time steps,               *
+!          2. Divergence of time-mean wind between the 2 time steps            *
+!*******************************************************************************
+! Variables:                                                                   *
+!                                                                              *
+! uuln, vvln    nested horiz. wind components on original eta levels           *
+!                                                                              *
+!                                                                              *
+!*******************************************************************************
+
+      use par_mod
+      use com_mod
+
+      integer itime
+
+      integer ix,jy,kz,k,n,ixp,jyp,ixm,jym,l
+      real pih,ylat,ylatp,ylatm
+      parameter(pih=pi/180.)
+      real xmasscolumn(2),pint(nzmax,2)
+      real pdiff(0:nxmaxn-1,0:nymaxn-1,nzmax,2)
+      real gradx(2),grady(2),gradterm
+      real divx,divy,div(nzmax,3),divwater(nzmax,3)
+      real divxwater,divywater,gradtermwater
+      real watercolumn(2),watermass(nzmax,2)
+      real dwatercolumndt(0:nxmaxn-1,0:nymaxn-1)
+      real divwatercolumn(0:nxmaxn-1,0:nymaxn-1,3)
+
+!      write(*,*) 'maxis: ',nxn(1),nyn(1),nxmaxn,nymaxn
+      do l=1,numbnests
+!     write(89,*) xlon0+dx,ylat0+dy,nxn(l)-2,nyn(l)-2,dx,dy,1,10000.,3
+
+! Loop over the whole grid
+!*************************
+
+      do 10 jy=0,nyn(l)-1
+        do 10 ix=0,nxn(l)-1
+
+!*****************************************************************
+! 1. Compute the precipitable water tendency in atmospheric column
+!    and the fluxes needed for the divergence term (step 2)
+!*****************************************************************
+
+          do 20 k=1,2
+            n=memind(k)
+
+            watercolumn(n)=0.
+
+! Compute mass of total atmospheric column (xmasscolumn)
+! and water vapor mass in atmospheric column (watercolumn)
+!*********************************************************
+
+            xmasscolumn(n)=psn(ix,jy,1,n,l)/ga
+
+            do 21 kz=1,nwz
+21            pint(kz,n)=akm(kz)+bkm(kz)*psn(ix,jy,1,n,l)
+
+            do 20 kz=2,nwz
+              pdiff(ix,jy,kz,n)=pint(kz-1,n)-pint(kz,n)
+              watermass(kz,n)=qvhn(ix,jy,kz,n,l)*pdiff(ix,jy,kz,n)/ga
+20            watercolumn(n)=watercolumn(n)+watermass(kz,n)
+
+
+! Compute the tendency
+!*********************
+
+          dwatercolumndt(ix,jy)=(watercolumn(memind(2))- &
+             watercolumn(memind(1)))/float(memtime(2)-memtime(1))
+
+
+10        continue
+
+
+!**************************************************
+! 2. Compute the precipitable water divergence term
+!**************************************************
+
+      do 50 jy=1,nyn(l)-2
+        ylat=ylat0n(l)+float(jy)*dyn(l)
+        jyp=jy+1
+        jym=jy-1
+        ylatp=ylat0n(l)+float(jyp)*dyn(l)
+        ylatm=ylat0n(l)+float(jym)*dyn(l)
+        do 50 ix=1,nxn(l)-1
+          ixp=ix+1
+          if (ixp.gt.nxn(l)-1) ixp=0
+          ixm=ix-1
+
+! Calculate gradient of surface pressure
+!***************************************
+
+          do 54 k=1,2
+            n=memind(k)
+            gradx(n)=(psn(ixp,jy,1,n,l)-psn(ixm,jy,1,n,l)) &
+              /(2.*dxn(l)*pih)/r_earth/cos(ylat*pih)
+            grady(n)=(psn(ix,jyp,1,n,l)-psn(ix,jym,1,n,l)) &
+               /(2.*dyn(l)*pih)/r_earth
+54          continue
+
+
+          do 51 kz=2,nwz
+            do 52 k=1,2
+              n=memind(k)
+
+              divx=(uuln(ixp,jy,kz,n,l)-uuln(ixm,jy,kz,n,l))/ &
+              (2.*dxn(l)*pih)/r_earth/cos(ylat*pih)
+              divy=(vvln(ix,jyp,kz,n,l)*cos(ylatp*pih)- &
+              vvln(ix,jym,kz,n,l)*cos(ylatm*pih))/ &
+              (2.*dyn(l)*pih)/r_earth/cos(ylat*pih)
+
+              gradterm=(uuln(ix,jy,kz,n,l)*gradx(n)+ &
+              vvln(ix,jy,kz,n,l)*grady(n))*(bkm(kz-1)-bkm(kz))
+              
+              div(kz,n)=(divx+divy)*pdiff(ix,jy,kz,n)+gradterm
+
+              divxwater=(uuln(ixp,jy,kz,n,l)*qvhn(ixp,jy,kz,n,l)- &
+                         uuln(ixm,jy,kz,n,l)*qvhn(ixm,jy,kz,n,l))/ &
+                        (2.*dxn(l)*pih)/r_earth/cos(ylat*pih)
+              divywater=(vvln(ix,jyp,kz,n,l)*qvhn(ix,jyp,kz,n,l)* &
+                     cos(ylatp*pih)- &
+                     vvln(ix,jym,kz,n,l)*qvhn(ix,jym,kz,n,l)* &
+                     cos(ylatm*pih))/&
+                     (2.*dyn(l)*pih)/r_earth/cos(ylat*pih)
+
+              gradtermwater=(uuln(ix,jy,kz,n,l)*qvhn(ix,jy,kz,n,l) &
+                *gradx(n)+ &
+                vvln(ix,jy,kz,n,l)*qvhn(ix,jy,kz,n,l)*grady(n))* &
+               (bkm(kz-1)-bkm(kz))
+              
+              divwater(kz,n)=(divxwater+divywater)*pdiff(ix,jy,kz,n)+ &
+                gradtermwater
+
+52            continue
+
+            div(kz,3)=(div(kz,1)+div(kz,2))/2.
+            divwater(kz,3)=(divwater(kz,1)+divwater(kz,2))/2.
+51          continue
+
+          do 56 n=1,3
+            divwatercolumn(ix,jy,n)=0.
+            do 53 kz=2,nwz
+53            divwatercolumn(ix,jy,n)=divwatercolumn(ix,jy,n)+ &
+                divwater(kz,n)/ga
+
+56          continue
+
+          n=memind(1)
+          e_minus_p_nest(ix,jy,n,l)=   &
+               dwatercolumndt(ix,jy)+divwatercolumn(ix,jy,3)
+
+! Convert to mm/day
+
+          e_minus_p_nest(ix,jy,n,l)=e_minus_p_nest(ix,jy,n,l)*86400.
+
+50        continue
+
+
+
+! Set values at the poles equal to values 1 deg equatorward
+!**********************************************************
+      n=memind(1)
+
+      do 80 ix=1,nxn(l)-1
+        e_minus_p_nest(ix,nyn(l)-1,n,l)=e_minus_p_nest(ix,nyn(l)-2,n,l)
+        e_minus_p_nest(ix,0,n,l)=e_minus_p_nest(ix,1,n,l)
+80      continue
+
+      write(90,'(i2,i8,1x,a11)') n,itime,' EminusP    '
+      do 85 jy=1,nyn(l)-2
+          write(90,*) (e_minus_p_nest(ix,jy,n,l),ix=1,nxn(l)-2)
+85        continue
+
+
+      end do
+      end