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