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