euler_rain.f90 5.77 KB
Newer Older
1
      subroutine euler_rain(itime)
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
!*******************************************************************************
! 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:                                                                   *
!                                                                              *
! uul, vvl    horiz. wind components on original eta levels                    *
!                                                                              *
!                                                                              *
!*******************************************************************************

      use par_mod
      use com_mod

17
18
      integer itime

19
20
21
22
23
24
25
26
27
28
      integer ix,jy,kz,k,n,ixp,jyp,ixm,jym
      real pih,ylat,ylatp,ylatm
      parameter(pih=pi/180.)
      real xmasscolumn(2),pint(nzmax,2)
      real pdiff(0:nxmax-1,0:nymax-1,nzmax,2),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:nxmax-1,0:nymax-1)
      real divwatercolumn(0:nxmax-1,0:nymax-1,3)
29
      integer nmax
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
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
153
154
155
156
157
158
159
160
161
162
163
      parameter(nmax=nxmax*nymax)

      write(89,*) xlon0+dx,ylat0+dy,nx-2,ny-2,dx,dy,1,10000.,3

! Loop over the whole grid
!*************************

      do 10 jy=0,ny-1
        do 10 ix=0,nx-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)=ps(ix,jy,1,n)/ga

            do 21 kz=1,nwz
21            pint(kz,n)=akm(kz)+bkm(kz)*ps(ix,jy,1,n)

            do 20 kz=2,nwz
              pdiff(ix,jy,kz,n)=pint(kz-1,n)-pint(kz,n)
              watermass(kz,n)=qvh(ix,jy,kz,n)*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,ny-2
        ylat=ylat0+float(jy)*dy
        jyp=jy+1
        jym=jy-1
        ylatp=ylat0+float(jyp)*dy
        ylatm=ylat0+float(jym)*dy
        do 50 ix=1,nx-1
          ixp=ix+1
          if (ixp.gt.nx-1) ixp=0
          ixm=ix-1

! Calculate gradient of surface pressure
!***************************************

          do 54 k=1,2
            n=memind(k)
            gradx(n)=(ps(ixp,jy,1,n)-ps(ixm,jy,1,n))/(2.*dx*pih)/ &
            r_earth/cos(ylat*pih)
            grady(n)=(ps(ix,jyp,1,n)-ps(ix,jym,1,n))/(2.*dy*pih)/ &
            r_earth
54          continue


          do 51 kz=2,nwz
            do 52 k=1,2
              n=memind(k)

              divx=(uul(ixp,jy,kz,n)-uul(ixm,jy,kz,n))/ &
              (2.*dx*pih)/r_earth/cos(ylat*pih)
              divy=(vvl(ix,jyp,kz,n)*cos(ylatp*pih)- &
              vvl(ix,jym,kz,n)*cos(ylatm*pih))/ &
              (2.*dy*pih)/r_earth/cos(ylat*pih)

              gradterm=(uul(ix,jy,kz,n)*gradx(n)+ &
              vvl(ix,jy,kz,n)*grady(n))*(bkm(kz-1)-bkm(kz))
              
              div(kz,n)=(divx+divy)*pdiff(ix,jy,kz,n)+gradterm

              divxwater=(uul(ixp,jy,kz,n)*qvh(ixp,jy,kz,n)- &
                         uul(ixm,jy,kz,n)*qvh(ixm,jy,kz,n))/ &
                        (2.*dx*pih)/r_earth/cos(ylat*pih)
              divywater=(vvl(ix,jyp,kz,n)*qvh(ix,jyp,kz,n)* &
                     cos(ylatp*pih)- &
                     vvl(ix,jym,kz,n)*qvh(ix,jym,kz,n)*cos(ylatm*pih))/&
                     (2.*dy*pih)/r_earth/cos(ylat*pih)

              gradtermwater=(uul(ix,jy,kz,n)*qvh(ix,jy,kz,n)*gradx(n)+ &
                vvl(ix,jy,kz,n)*qvh(ix,jy,kz,n)*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

          e_minus_p(ix,jy)=dwatercolumndt(ix,jy)+divwatercolumn(ix,jy,3)

! Convert to mm/day

          e_minus_p(ix,jy)=e_minus_p(ix,jy)*86400.

50        continue



! Set values at the poles equal to values 1 deg equatorward
!**********************************************************

      do 80 ix=1,nx-1
        e_minus_p(ix,ny-1)=e_minus_p(ix,ny-2)
        e_minus_p(ix,0)=e_minus_p(ix,1)
80      continue

164
      write(89,'(i8,1x,a11)') itime,' EminusP    '
165
      do 85 jy=1,ny-2
166
          write(89,*) (e_minus_p(ix,jy),ix=1,nx-2)
167
168
169
170
171
172
173
174
175
176
177
85        continue

! Use results only if we are inside the simulation period (don't use the first fields,
! which are probably only needed for interpolation to time 0
!*************************************************************************************

      if (memtime(2).gt.0) then
      endif


      end