UPP (develop)
Loading...
Searching...
No Matches
CALTAU.f
Go to the documentation of this file.
1
27!-----------------------------------------------------------------------
32!-----------------------------------------------------------------------
33 SUBROUTINE caltau(TAUX,TAUY)
34
35!
36!
37 use vrbls3d, only: zint, pmid, q, t, uh, vh, el_pbl, zmid
38 use vrbls2d, only: z0, uz0, vz0
39 use masks, only: lmh
40 use params_mod, only: d00, d50, h1, d608, rd, d25
41 use ctlblk_mod, only: jsta_2l, jend_2u, lm, jsta, jend, spval, jsta_m,&
42 jm, im, jend_m, ista, iend, ista_m, iend_m, ista_2l, iend_2u
43 use gridspec_mod, only: gridtype
44
45 implicit none
46!
47! DECLARE VARIABLES.
48 INTEGER, dimension(4) :: KK(4)
49 INTEGER, dimension(jm) :: ive, ivw
50 REAL, dimension(ista:iend,jsta:jend), intent(inout) :: TAUX, TAUY
51 REAL, ALLOCATABLE :: EL(:,:,:)
52 REAL, dimension(ista:iend,jsta:jend) :: EGRIDU,EGRIDV,EGRID4,EGRID5, EL0
53 REAL UZ0V,VZ0V
54 CHARACTER*1 AGRID
55 integer I,J,LMHK,IE,IW,ii,jj
56 real DZ,RDZ,RSFC,TV,RHO,ULMH,VLMH,DELUDZ,DELVDZ,ELSQR,ZINT1, &
57 zint2,z0v,psfc,tvv,qvv,elv,elv1,elv2
58!
59!********************************************************************
60! START CALTAU HERE.
61!
62 ALLOCATE (el(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
63!
64! COMPUTE MASTER LENGTH SCALE.
65!
66! CALL CLMAX(EL0,EGRIDU,EGRIDV,EGRID4,EGRID5)
67! CALL MIXLEN(EL0,EL)
68!
69! INITIALIZE OUTPUT AND WORK ARRAY TO ZERO.
70!
71 DO j=jsta,jend
72 DO i=ista,iend
73 egridu(i,j) = d00
74 egridv(i,j) = d00
75 taux(i,j) = spval
76 tauy(i,j) = spval
77 ENDDO
78 ENDDO
79!
80! COMPUTE SURFACE LAYER U AND V WIND STRESSES.
81!
82! ASSUME THAT U AND V HAVE UPDATED HALOS
83!
84 IF(gridtype == 'A')THEN
85 CALL clmax(el0,egridu,egridv,egrid4,egrid5)
86 CALL mixlen(el0,el)
87
88 DO j=jsta,jend
89 DO i=ista,iend
90!
91 lmhk = nint(lmh(i,j))
92 IF(el(i,j,lmhk-1)<spval.and.z0(i,j)<spval.and. &
93 uz0(i,j)<spval.and.vz0(i,j)<spval)THEN
94!
95! COMPUTE THICKNESS OF LAYER AT MASS POINT.
96!
97 dz = d50*(zint(i,j,lmhk)-zint(i,j,lmhk+1))
98 dz = dz-z0(i,j)
99 rdz = 1./dz
100!
101! COMPUTE REPRESENTATIVE AIR DENSITY.
102!
103 psfc = pmid(i,j,lmhk)
104 tv = (h1+d608*q(i,j,lmhk))*t(i,j,lmhk)
105 rho = psfc/(rd*tv)
106!
107! COMPUTE A MEAN MASS POINT WIND IN THE
108! FIRST ATMOSPHERIC ETA LAYER.
109!
110 ulmh = uh(i,j,lmhk)
111 vlmh = vh(i,j,lmhk)
112!
113! COMPUTE WIND SHEAR COMPONENTS ACROSS LAYER.
114!
115 deludz = (ulmh-uz0(i,j))*rdz
116 delvdz = (vlmh-vz0(i,j))*rdz
117!
118! COMPUTE U (EGRIDU) AND V (EGRIDV) WIND STRESSES.
119!
120 elsqr = el(i,j,lmhk-1)*el(i,j,lmhk-1)
121 taux(i,j) = rho*elsqr*deludz*deludz
122 tauy(i,j) = rho*elsqr*delvdz*delvdz
123 ELSE
124 taux(i,j) = spval
125 tauy(i,j) = spval
126 ENDIF
127
128!
129 END DO
130 END DO
131 ELSE IF(gridtype == 'E')THEN
132 call exch(zint(1,jsta_2l,lm))
133 call exch(zint(1,jsta_2l,lm+1))
134 call exch(z0(1,jsta_2l))
135 call exch(pmid(1,jsta_2l,lm))
136 call exch(t(1,jsta_2l,lm))
137 call exch(q(1,jsta_2l,lm))
138 call exch(el_pbl(1,jsta_2l,lm))
139 call exch(el_pbl(1,jsta_2l,lm-1))
140
141 DO j=jsta_m,jend_m
142 ive(j)=mod(j,2)
143 ivw(j)=ive(j)-1
144 ENDDO
145
146 DO j=jsta_m,jend_m
147 DO i=ista_m,iend_m
148!
149 lmhk = nint(lmh(i,j))
150 ie=i+ive(j)
151 iw=i+ivw(j)
152 zint1=(zint(iw,j,lmhk)+zint(ie,j,lmhk) &
153 +zint(i,j+1,lmhk)+zint(i,j-1,lmhk))*d25
154 zint2=(zint(iw,j,lmhk+1)+zint(ie,j,lmhk+1) &
155 +zint(i,j+1,lmhk+1)+zint(i,j-1,lmhk+1))*d25
156 dz = d50*(zint1-zint2)
157 z0v=(z0(iw,j)+z0(ie,j)+z0(i,j+1)+z0(i,j-1))*d25
158 dz = dz-z0v
159 rdz = 1./dz
160!
161! COMPUTE REPRESENTATIVE AIR DENSITY.
162!
163 psfc = (pmid(iw,j,lmhk)+pmid(ie,j,lmhk) &
164 +pmid(i,j+1,lmhk)+pmid(i,j-1,lmhk))*d25
165 tvv = (t(iw,j,lmhk)+t(ie,j,lmhk) &
166 +t(i,j+1,lmhk)+t(i,j-1,lmhk))*d25
167 qvv = (q(iw,j,lmhk)+q(ie,j,lmhk) &
168 +q(i,j+1,lmhk)+q(i,j-1,lmhk))*d25
169 tv = (h1+d608*qvv)*tvv
170 rho = psfc/(rd*tv)
171
172! COMPUTE WIND SHEAR COMPONENTS ACROSS LAYER.
173!
174 deludz = (uh(i,j,lmhk)-uz0(i,j))*rdz
175 delvdz = (vh(i,j,lmhk)-vz0(i,j))*rdz
176
177! COMPUTE U (EGRIDU) AND V (EGRIDV) WIND STRESSES.
178!
179 elv1=(el_pbl(iw,j,lmhk)+el_pbl(ie,j,lmhk) &
180 +el_pbl(i,j+1,lmhk)+el_pbl(i,j-1,lmhk))*d25
181 elv2=(el_pbl(iw,j,lmhk-1)+el_pbl(ie,j,lmhk-1) &
182 +el_pbl(i,j+1,lmhk-1)+el_pbl(i,j-1,lmhk-1))*d25
183 elv=(elv1+elv2)/2.0 ! EL is defined at the bottom of layer
184 elsqr =elv*elv
185 taux(i,j)=rho*elsqr*deludz*deludz
186 tauy(i,j)=rho*elsqr*delvdz*delvdz
187! ii=im/2
188! jj=(jsta+jend)/2
189! if(i==ii.and.j==jj)print*,'sample tau'
190! & ,RHO,ELSQR,DELUDZ,DELVDZ
191 END DO
192 END DO
193 ELSE IF(gridtype == 'B')THEN
194! PUT TAUX AND TAUY ON MASS POINTS
195 call exch(vh(1,jsta_2l,lm))
196 DO j=jsta_m,jend_m
197 DO i=ista_m,iend_m
198!
199 lmhk = nint(lmh(i,j))
200!
201! COMPUTE THICKNESS OF LAYER AT MASS POINT.
202!
203! DZ = D50*(ZINT(I,J,LMHK)-ZINT(I,J,LMHK+1))
204! DZ = ZMID(I,J,LMHK)-Z0(I,J)
205 dz=zmid(i,j,lmhk)-(z0(i,j)+zint(i,j,lmhk+1))
206 if(dz==0.0)dz=0.2
207 rdz = 1./dz
208!
209! COMPUTE REPRESENTATIVE AIR DENSITY.
210!
211 psfc = pmid(i,j,lmhk)
212 tv = (h1+d608*q(i,j,lmhk))*t(i,j,lmhk)
213 rho = psfc/(rd*tv)
214!
215! PUT U AND V ONTO MASS POINTS
216!
217 ulmh = 0.5*(uh(i-1,j,lmhk)+uh(i,j,lmhk))
218 vlmh = 0.5*(vh(i,j-1,lmhk)+vh(i,j,lmhk))
219!
220! COMPUTE WIND SHEAR COMPONENTS ACROSS LAYER.
221!
222 deludz = (ulmh-uz0(i,j))*rdz
223 delvdz = (vlmh-vz0(i,j))*rdz
224!
225! COMPUTE U (EGRIDU) AND V (EGRIDV) WIND STRESSES.
226!
227 elv=0.5*(el_pbl(i,j,lmhk)+el_pbl(i,j,lmhk-1))
228 elsqr = elv*elv
229 taux(i,j) = rho*elsqr*deludz*deludz
230! if(TAUX(I,J)>1.0e2)print*,'Debug TAUX= ',i,j, &
231! ELV,ULMH,UZ0(I,J),ZMID(I,J,LMHK),Z0(I,J),RDZ,TAUX(I,J),zint(i,j,lm+1)
232 tauy(i,j) = rho*elsqr*delvdz*delvdz
233
234 END DO
235 END DO
236 END IF
237!
238 DEALLOCATE(el)
239! END OF ROUTINE.
240!
241 RETURN
242 END
subroutine caltau(taux, tauy)
Subroutine that computes U and V wind stresses.
Definition CALTAU.f:34
subroutine exch(a)
exch() Subroutine that exchanges one halo row.
Definition EXCH.f:27