euler_rain.f90 5.52 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
29
      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)

Sabine's avatar
Sabine committed
30
!     write(89,*) xlon0+dx,ylat0+dy,nx-2,ny-2,dx,dy,1,10000.,3
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

! 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

Sabine's avatar
Sabine committed
144
145
146
          n=memind(1)
          e_minus_p(ix,jy,n)=   &
               dwatercolumndt(ix,jy)+divwatercolumn(ix,jy,3)
147
148
149

! Convert to mm/day

Sabine's avatar
Sabine committed
150
          e_minus_p(ix,jy,n)=e_minus_p(ix,jy,n)*86400.
151
152
153
154
155
156
157

50        continue



! Set values at the poles equal to values 1 deg equatorward
!**********************************************************
Sabine's avatar
Sabine committed
158
      n=memind(1)
159
160

      do 80 ix=1,nx-1
Sabine's avatar
Sabine committed
161
162
        e_minus_p(ix,ny-1,n)=e_minus_p(ix,ny-2,n)
        e_minus_p(ix,0,n)=e_minus_p(ix,1,n)
163
164
80      continue

Sabine's avatar
Sabine committed
165
      write(89,'(i2,i8,1x,a11)') n,itime,' EminusP    '
166
      do 85 jy=1,ny-2
Sabine's avatar
Sabine committed
167
          write(89,*) (e_minus_p(ix,jy,n),ix=1,nx-2)
168
169
170
85        continue

      end