CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
module_SF_JSFC.F90
1!-----------------------------------------------------------------------
2!
4!
5!-----------------------------------------------------------------------
6!
7!*** THE J SURFACE SCHEME
8!
9!-----------------------------------------------------------------------
10!
11! USE MODULE_INCLUDE
12!
13! USE MODULE_CONSTANTS,ONLY : A2,A3,A4,CP,ELWV &
14! ,G,P608,PI,PQ0,R_D,R_V,CAPPA
15!
16!-----------------------------------------------------------------------
17!
18
19 USE machine, only: kfpt => kind_phys
20
21 IMPLICIT NONE
22!
23!-----------------------------------------------------------------------
24!
25! integer,parameter :: isingle=selected_int_kind(r=9)
26! integer,parameter :: idouble=selected_int_kind(r=18)
27! integer,parameter :: single=selected_real_kind(p=6,r=37)
28! integer,parameter :: double=selected_real_kind(p=13,r=200)
29!
30! integer,parameter:: &
31! klog=4 &
32! ,kint=isingle &
33! ,kdin=idouble &
34! ,kfpt=single &
35! ,kdbl=double
36
37! real (kind=kfpt),parameter :: r4_in=x'ffbfffff'
38! real (kind=kdbl),parameter :: r8_in=x'fff7ffffffffffff'
39! integer(kind=kint),parameter :: i4_in=-999 ! -huge(1)
40!
41 ! integer,parameter:: &
42 ! klog=4 & ! logical variables
43 ! ,kint=4 & ! integer variables
44 ! !,kfpt=4 & ! floating point variables
45 ! ,kfpt=8 & ! floating point variables
46 ! ,kdbl=8 ! double precision
47!
48 PRIVATE
49!
50 PUBLIC :: jsfc_init,jsfc
51!
52 INTEGER :: itrmx=5 ! Iteration count for mixing length computation
53!
54 REAL(kind=kfpt),PARAMETER :: vkarman=0.4
55
56 REAL(kind=kfpt),PARAMETER :: a2=17.2693882,a3=273.15,a4=35.86,cp=1004.6 &
57 ,elwv=2.501e6,epsq2=0.02,g=9.8060226 &
58 ,pq0=379.90516,r_d=287.04,r_v=461.6 &
59 ,p608=r_v/r_d-1.,cappa=r_d/cp &
60 ,pi=3.141592653589793
61
62 REAL(kind=kfpt),PARAMETER :: xlv=elwv
63 REAL(kind=kfpt),PARAMETER :: elocp=2.72e6/cp
64 REAL(kind=kfpt),PARAMETER :: a2s=17.2693882,a3s=273.16,a4s=35.86
65 REAL(kind=kfpt),PARAMETER :: glkbr=10.,glkbs=30. &
66 ,qvisc=2.1e-5,ric=0.505,small=0.35 &
67 ,sqpr=0.84,sqsc=0.84,sqvisc=258.2 &
68 ,tvisc=2.1e-5 &
69 ,ustc=0.7,ustr=0.225,visc=1.5e-5 &
70 ,wwst=1.2,ztfc=1.
71 REAL(kind=kfpt),PARAMETER :: seafc=0.98,pq0sea=pq0*seafc
72
73 REAL(kind=kfpt),PARAMETER :: cziv=small*glkbs,grrs=glkbr/glkbs
74
75 REAL(kind=kfpt),PARAMETER :: rtvisc=1./tvisc,rvisc=1./visc &
76 ,zqrzt=sqsc/sqpr
77
78 REAL(kind=kfpt),PARAMETER :: ustfc=0.018/g &
79 ,fzq1=rtvisc*qvisc*zqrzt &
80 ,fzq2=rtvisc*qvisc*zqrzt &
81 ,fzt1=rvisc *tvisc*sqpr &
82 ,fzt2=cziv*grrs*tvisc*sqpr &
83 ,fzu1=cziv*visc
84 REAL(kind=kfpt),PARAMETER :: wwst2=wwst*wwst &
85 ,rqvisc=1./qvisc
86
87 REAL(kind=kfpt),PARAMETER :: rcap=1./cappa
88 REAL(kind=kfpt),PARAMETER :: gocp02=g/cp*2.,gocp10=g/cp*10.
89 REAL(kind=kfpt),PARAMETER :: epsu2=1.e-6,epsust=1.e-9,epszt=1.e-28
90 REAL(kind=kfpt),PARAMETER :: czil=0.1,excml=0.0001,excms=0.0001 &
91 & ,fh=1.10,topofac=9.0e-6
92
93 REAL(kind=kfpt),PARAMETER :: zilfc=-czil*vkarman*sqvisc
94 REAL(kind=kfpt),PARAMETER :: epsq=1.e-9
95!
96!-----------------------------------------------------------------------
97 INTEGER, PARAMETER :: kztm=10001,kztm2=kztm-2
98!
99 REAL(kind=kfpt),PRIVATE,SAVE :: &
100 dzeta1,dzeta2,fh01,fh02,ztmax1,ztmax2,ztmin1,ztmin2
101!
102 REAL(kind=kfpt),DIMENSION(KZTM),PRIVATE,SAVE :: &
103 psih1,psih2,psim1,psim2
104!
105 INTEGER :: ierr
106!
107!-----------------------------------------------------------------------
108!
109 CONTAINS
110!
111!-----------------------------------------------------------------------
112!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
113!-----------------------------------------------------------------------
114 SUBROUTINE jsfc(FLAG_ITER,ITER,ME &
115 & ,NTSD,EPSQ2,HT,DZ &
116 & ,PHMID,PHINT,TH,T,Q,QC,U,V,Q2 &
117 & ,TSK,QSFC,THZ0,QZ0,UZ0,VZ0 &
118 & ,XLAND &
119 & ,USTAR,Z0,Z0BASE,PBLH,MAVAIL,RMOL &
120 & ,AKHS,AKMS,CHKLOWQ,HLFLX,RIB &
121 & ,CM,CH,STRESS,FFM,FFH,WIND,FM10,FH2 &
122 & ,A1U,A1T,A1Q &
123 & ,IDS,IDE,JDS,JDE,KDS,KDE &
124 & ,IMS,IME,JMS,JME,KMS,KME &
125 & ,ITS,ITE,JTS,JTE,KTS,LM,errmsg,errflg)
126!
127!-----------------------------------------------------------------------
128! SUBROUTINE JSFC(NTSD,EPSQ2,HT,DZ &
129! & ,PHMID,PHINT,TH,T,Q,QC,U,V,Q2 &
130! & ,TSK,QSFC,THZ0,QZ0,UZ0,VZ0 &
131! & ,XLAND &
132! & ,VEGFRC,SNOWC & !added 5/17/2013
133! & ,USTAR,Z0,Z0BASE,PBLH,MAVAIL,RMOL &
134! & ,AKHS,AKMS &
135! & ,CHS,CHS2,CQS2,HFX,QFX,FLX_LH,FLHC,FLQC &
136! & ,QGH,CPM,CT &
137! & ,U10,V10,T02,TH02,TSHLTR,TH10,Q02,QSHLTR,Q10 &
138! & ,PSHLTR,RIB & ! Added Bulk Richardson No.
139! & ,IDS,IDE,JDS,JDE,KDS,KDE &
140! & ,IMS,IME,JMS,JME,KMS,KME &
141! & ,ITS,ITE,JTS,JTE,KTS,LM)
142!----------------------------------------------------------------------
143!
144 IMPLICIT NONE
145!
146!----------------------------------------------------------------------
147 INTEGER,INTENT(IN) :: ids,ide,jds,jde,kds,kde &
148 & ,IMS,IME,JMS,JME,KMS,KME &
149 & ,ITS,ITE,JTS,JTE,KTS,LM
150!
151 INTEGER,INTENT(IN) :: ntsd,iter,me
152 LOGICAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: flag_iter
153 real(kind=kfpt),dimension(1:lm),intent(in):: epsq2
154!
155 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: ht,mavail,tsk &
156 & ,xland,z0base
157! & ,VEGFRC,SNOWC
158!
159 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(IN) :: dz,phmid
160!
161 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME,1:LM+1),INTENT(IN) :: phint
162!
163 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(IN) :: q,qc,u,v,q2,t,th
164!
165! REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: FLX_LH,HFX,PSHLTR &
166! & ,QFX,Q10,QSHLTR &
167! & ,TH10,TSHLTR,T02 &
168! & ,U10,V10,TH02,Q02
169 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME) :: flx_lh,hfx &
170 & ,qfx,q10,th10,t02 &
171 & ,u10,v10,th02,q02
172 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME) :: pshltr,qshltr,tshltr
173!
174 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: akhs,akms &
175 & ,pblh,qsfc,rib
176!
177 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: qz0,rmol,thz0 &
178 & ,ustar,uz0,vz0 &
179 & ,z0
180!
181 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: hlflx,chklowq
182 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: cm,ch,stress,ffm &
183 & ,ffh,wind,fm10,fh2 &
184 & ,a1u,a1t,a1q
185 character(len=*),intent(out) :: errmsg
186 integer, intent(out) :: errflg
187!
188! REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: CHS,CHS2,CQS2 &
189! & ,CPM,CT,FLHC,FLQC &
190! & ,QGH
191 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME) :: chs,chs2,cqs2 &
192 & ,flhc,flqc
193 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME) :: qgh,cpm,ct
194!----------------------------------------------------------------------
195!***
196!*** LOCAL VARIABLES
197!***
198 INTEGER :: i,j,k,lmh,lpbl
199!
200 REAL(kind=kfpt) :: a,apesfc,b,btgx,cwmlow &
201 & ,dqdt,dtdif,dtdt,dudt,dvdt &
202 & ,fis &
203 & ,p02p,p10p,plow,psfc,ptop,qlow,qs02,qs10 &
204 & ,rapa,rapa02,rapa10,ratiomx,rdz,seamask,sm &
205 & ,t02p,t10p,tem,th02p,th10p,thlow,thelow,thm &
206 & ,tlow,tz0,ulow,vlow,zsl
207!
208 REAL(kind=kfpt),DIMENSION(KTS:LM) :: cwmk,pk,q2k,qk,thek,thk,tk,uk,vk
209!
210 REAL(kind=kfpt),DIMENSION(KTS:LM-1) :: el,elm
211!
212 REAL(kind=kfpt),DIMENSION(KTS:LM+1) :: zhk
213!
214 REAL(kind=kfpt),DIMENSION(ITS:ITE,JTS:JTE) :: thsk
215!
216 REAL(kind=kfpt),DIMENSION(ITS:ITE,JTS:JTE,KTS:LM+1) :: zint
217!
218!----------------------------------------------------------------------
219!**********************************************************************
220 ! Initialize error-handling
221 errflg = 0
222 errmsg = ''
223!----------------------------------------------------------------------
224!
225!*** MAKE PREPARATIONS
226!
227!----------------------------------------------------------------------
228 DO j=jts,jte
229 DO i=its,ite
230 IF(flag_iter(i,j))THEN
231 DO k=kts,lm+1
232 zint(i,j,k)=0.
233 ENDDO
234 END IF
235 ENDDO
236 ENDDO
237!
238 DO j=jts,jte
239 DO i=its,ite
240 IF(flag_iter(i,j))THEN
241 zint(i,j,lm+1)=ht(i,j) ! Z at bottom of lowest sigma layer
242 pblh(i,j)=-1.
243!
244!!!!!!!!!
245!!!!!! UNCOMMENT THESE LINES IF USING ETA COORDINATES
246!!!!!!!!!
247!!!!!! ZINT(I,J,LM+1)=1.E-4 ! Z of bottom of lowest eta layer
248!!!!!! ZHK(LM+1)=1.E-4 ! Z of bottom of lowest eta layer
249!
250 END IF
251 ENDDO
252 ENDDO
253!
254 DO j=jts,jte
255 DO i=its,ite
256 IF(flag_iter(i,j))THEN
257 DO k=lm,kts,-1
258 zint(i,j,k)=zint(i,j,k+1)+dz(i,j,k)
259 ENDDO
260 END IF
261 ENDDO
262 ENDDO
263!
264 IF(ntsd==0.and.iter==1) then
265 DO j=jts,jte
266 DO i=its,ite
267 IF(flag_iter(i,j))THEN
268 ustar(i,j)=0.1
269 fis=ht(i,j)*g
270 sm=xland(i,j)-1.
271!!! Z0(I,J)=SM*Z0SEA+(1.-SM)*(Z0(I,J)*Z0MAX+FIS*FCM+Z0LAND)
272!!! Z0(I,J)=SM*Z0SEA+(1.-SM)*(Z0(I,J)*Z0MAX+FIS*FCM+Z0LAND)
273 END IF
274 ENDDO
275 ENDDO
276 ENDIF
277!
278!!!! IF(NTSD==1)THEN
279 DO j=jts,jte
280 DO i=its,ite
281 ct(i,j)=0.
282 ENDDO
283 ENDDO
284!!!! ENDIF
285!
286!......................................................................
287!$omp parallel do &
288!$omp private (j,i,lmh,ptop,psfc,seamask,k,thk,tk,ratiomx,qk,pk, &
289!$omp cwmk,thek,q2k,zhk,uk,vk,lpbl,plow,tlow,thlow,thelow, &
290!$omp qlow,cwmlow,ulow,vlow,zsl,apesfc,tz0,rapa,th02p,th10p, &
291!$omp rapa02,rapa10,t02p,t10p,p02p,p10p,qs02,qs10)
292!......................................................................
293!----------------------------------------------------------------------
294 setup_integration: DO j=jts,jte
295!----------------------------------------------------------------------
296!
297 DO i=its,ite
298 IF(flag_iter(i,j))THEN
299!
300!*** LOWEST LAYER ABOVE GROUND MUST BE FLIPPED
301!
302 lmh=lm
303!
304 ptop=phint(i,j,1)
305 psfc=phint(i,j,lmh+1)
306! Define THSK here (for first timestep mostly)
307 thsk(i,j)=tsk(i,j)/(psfc*1.e-5)**cappa
308!
309!*** CONVERT LAND MASK (1 FOR SEA; 0 FOR LAND)
310!
311 seamask=xland(i,j)-1.
312!
313!*** FILL 1-D VERTICAL ARRAYS
314!
315 DO k=lm,kts,-1
316 thk(k)=th(i,j,k)
317 tk(k)=t(i,j,k)
318 qk(k)=q(i,j,k)
319 pk(k)=phmid(i,j,k)
320 cwmk(k)=qc(i,j,k)
321 thek(k)=(cwmk(k)*(-elocp/tk(k))+1.)*thk(k)
322 q2k(k)=q2(i,j,k)
323!
324!
325!*** COMPUTE THE HEIGHTS OF THE LAYER INTERFACES
326!
327 zhk(k)=zint(i,j,k)
328!
329 ENDDO
330 zhk(lm+1)=ht(i,j) ! Z at bottom of lowest sigma layer
331!
332 DO k=lm,kts,-1
333 uk(k)=u(i,j,k)
334 vk(k)=v(i,j,k)
335 ENDDO
336!
337!*** FIND THE HEIGHT OF THE PBL
338!
339 lpbl=lmh
340 DO k=lmh-1,1,-1
341 IF(q2k(k)<=epsq2(k)*fh) THEN
342 lpbl=k
343 GO TO 110
344 ENDIF
345 ENDDO
346!
347 lpbl=1
348!
349!-----------------------------------------------------------------------
350!--------------THE HEIGHT OF THE PBL------------------------------------
351!-----------------------------------------------------------------------
352!
353 110 pblh(i,j)=zhk(lpbl)-zhk(lmh+1)
354!
355!----------------------------------------------------------------------
356 IF(qc(i,j,lm).GT.epsq)THEN
357 chklowq(i,j)=0.
358 ELSE
359 chklowq(i,j)=1.
360 ENDIF
361!***
362!*** FIND THE SURFACE EXCHANGE COEFFICIENTS
363!***
364!----------------------------------------------------------------------
365 plow=pk(lmh)
366 tlow=tk(lmh)
367 thlow=thk(lmh)
368 thelow=thek(lmh)
369 qlow=qk(lmh)
370 cwmlow=cwmk(lmh)
371 ulow=uk(lmh)
372 vlow=vk(lmh)
373 zsl=(zhk(lmh)-zhk(lmh+1))*0.5
374! if(me.eq.0)print*,'ZSL,ZHK(LMH),ZHK(LMH+1,LMH=',ZSL,ZHK(LMH),ZHK(LMH+1),LMH
375 apesfc=(psfc*1.e-5)**cappa
376 if(ntsd==0) then
377 tz0=tsk(i,j)
378 else
379 tz0=thz0(i,j)*apesfc
380 endif
381!
382 CALL sfcdif(ntsd,seamask,thsk(i,j),qsfc(i,j),psfc &
383 & ,uz0(i,j),vz0(i,j),tz0,thz0(i,j),qz0(i,j) &
384 & ,ustar(i,j),z0(i,j),z0base(i,j),ct(i,j),rmol(i,j) &
385 & ,akms(i,j),akhs(i,j),pblh(i,j),mavail(i,j) &
386 & ,chs(i,j),chs2(i,j),cqs2(i,j) &
387 & ,hfx(i,j),qfx(i,j),flx_lh(i,j) &
388 & ,flhc(i,j),flqc(i,j),qgh(i,j),cpm(i,j) &
389 & ,ulow,vlow,tlow,thlow,thelow,qlow,cwmlow &
390 & ,zsl,plow,hlflx(i,j) &
391! & ,VEGFRC(I,J),SNOWC(I,J) & !added 5/17/2013
392 & ,u10(i,j),v10(i,j),tshltr(i,j),th10(i,j) &
393 & ,qshltr(i,j),q10(i,j),pshltr(i,j) &
394 & ,ffm(i,j),ffh(i,j),fm10(i,j),fh2(i,j) &
395 & ,a1u(i,j),a1t(i,j),a1q(i,j) &
396 & ,ids,ide,jds,jde,kds,kde &
397 & ,ims,ime,jms,jme,kms,kme &
398 & ,its,ite,jts,jte,kts,lm,i,j,zhk(lmh+1),rib(i,j) & ! Added Bulk Richardson No.
399 & ,errmsg, errflg)
400!
401!*** REMOVE SUPERATURATION AT 2M AND 10M
402!
403 rapa=apesfc
404 th02p=tshltr(i,j)
405 th10p=th10(i,j)
406 th02(i,j)=tshltr(i,j)
407!
408 rapa02=rapa-gocp02/th02p
409 rapa10=rapa-gocp10/th10p
410!
411 t02p=th02p*rapa02
412 t10p=th10p*rapa10
413! 1 may 06 tgs T02(I,J) = T02P
414 t02(i,j) = th02(i,j)*apesfc
415!
416 p02p=(rapa02**rcap)*1.e5
417 p10p=(rapa10**rcap)*1.e5
418!
419 qs02=pq0/p02p*exp(a2*(t02p-a3)/(t02p-a4))
420 qs10=pq0/p10p*exp(a2*(t10p-a3)/(t10p-a4))
421!
422 IF(qshltr(i,j)>qs02)qshltr(i,j)=qs02
423 IF(q10(i,j)>qs10)q10(i,j)=qs10
424 q02(i,j)=qshltr(i,j)/(1.-qshltr(i,j))
425!----------------------------------------------------------------------
426! STRESS(I,J)=USTAR(I,J)*USTAR(I,J)
427 wind(i,j)=max(ustar(i,j)*ffm(i,j)/vkarman,1.0)
428 cm(i,j)=vkarman*vkarman/(ffm(i,j)*ffm(i,j))
429 ch(i,j)=vkarman*vkarman/(ffm(i,j)*ffh(i,j))
430 tem=0.00001/dz(i,j,lm)
431 cm(i,j)=max(cm(i,j),tem)
432 ch(i,j)=max(ch(i,j),tem)
433 stress(i,j)=cm(i,j) * wind(i,j) * wind(i,j)
434 ustar(i,j)=sqrt(stress(i,j))
435!
436 END IF ! FLAG_ITER
437!
438 ENDDO
439!
440!----------------------------------------------------------------------
441!
442 ENDDO setup_integration
443!
444!......................................................................
445!$omp end parallel do
446!......................................................................
447!----------------------------------------------------------------------
448
449 END SUBROUTINE jsfc
450!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
451!----------------------------------------------------------------------
452 SUBROUTINE sfcdif(NTSD,SEAMASK,THS,QS,PSFC &
453 & ,UZ0,VZ0,TZ0,THZ0,QZ0 &
454 & ,USTAR,Z0,Z0BASE,CT,RLMO,AKMS,AKHS,PBLH,WETM &
455 & ,CHS,CHS2,CQS2,HFX,QFX,FLX_LH,FLHC,FLQC,QGH,CPM &
456 & ,ULOW,VLOW,TLOW,THLOW,THELOW,QLOW,CWMLOW &
457 & ,ZSL,PLOW,HLFLX &
458! & ,VEGF,SNOC & !added 5/17/2013
459 & ,U10,V10,TH02,TH10,Q02,Q10,PSHLTR &
460 & ,FFM,FFH,FM10,FH2,A1U,A1T,A1Q &
461 & ,IDS,IDE,JDS,JDE,KDS,KDE &
462 & ,IMS,IME,JMS,JME,KMS,KME &
463 & ,ITS,ITE,JTS,JTE,KTS,LM,I,J,ZSFC,RIB & ! Added Bulk Richardson No.
464 & ,errmsg, errflg)
465! ****************************************************************
466! * *
467! * SURFACE LAYER *
468! * *
469! ****************************************************************
470!----------------------------------------------------------------------
471!
472 IMPLICIT NONE
473!
474!----------------------------------------------------------------------
475 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
476 & ,IMS,IME,JMS,JME,KMS,KME &
477 & ,ITS,ITE,JTS,JTE,KTS,LM,i,j
478!
479 INTEGER,INTENT(IN) :: NTSD
480!
481 REAL(kind=kfpt),INTENT(IN) :: cwmlow,pblh,plow,qlow,psfc,seamask,zsfc &
482 & ,thelow,thlow,ths,tlow,tz0,ulow,vlow,wetm,zsl &
483 & ,z0base
484! ,VEGF,SNOC
485!
486 REAL(kind=kfpt),INTENT(OUT) :: chs,chs2,cpm,cqs2,ct,flhc,flqc,flx_lh,hfx &
487 & ,rib,pshltr,q02,q10,qfx,qgh,rlmo,th02,th10,u10,v10
488 REAL(kind=kfpt),INTENT(OUT) :: ffm,ffh,fm10,fh2,a1u,a1t,a1q
489!
490 REAL(kind=kfpt),INTENT(INOUT) :: akhs,akms,qz0,thz0,ustar,uz0,vz0,z0,qs
491 character(len=*), intent(out) :: errmsg
492 integer, intent(out) :: errflg
493!----------------------------------------------------------------------
494!***
495!*** LOCAL VARIABLES
496!***
497 INTEGER :: ITR,K
498!
499 REAL(kind=kfpt) :: a,b,btgh,btgx,cxchl,cxchs,dthv,du2,elfc,fct &
500 & ,hlflx,hsflx,hv,psh02,psh10,pshz,pshzl,psm10,psmz,psmzl &
501 & ,rdz,rdzt,rlma,rlmn,rlmp &
502 & ,rlogt,rlogu,rwgh,rz,rzst,rzsu,simh,simm,tem,thm &
503 & ,umflx,ustark,vmflx,wght,wghtt,wghtq,wstar2 &
504 & ,x,xlt,xlt4,xlu,xlu4,xt,xt4,xu,xu4,zetalt,zetalu &
505 & ,zetat,zetau,zq,zslt,zslu,zt,zu,topoterm,zzil,czetmax
506!vcw
507!
508!*** DIAGNOSTICS
509!
510 REAL(kind=kfpt) :: akhs02,akhs10,akms02,akms10,ekms10,qsat10,qsat2 &
511 & ,rlnt02,rlnt10,rlnu10,simh02,simh10,simm10,t02,t10 &
512 & ,term1,rlow,u10e,v10e,wstar,xlt02,xlt024,xlt10 &
513 & ,xlt104,xlu10,xlu104,xu10,xu104,zt02,zt10,ztat02,ztat10 &
514 & ,ztau,ztau10,zu10,zuuz
515! REAL(kind=kfpt) :: ZILFC1,SNOWZO, Zom_ztmax
516!----------------------------------------------------------------------
517!**********************************************************************
518!----------------------------------------------------------------------
519 ! Initialize error-handling
520 errflg = 0
521 errmsg = ''
522
523 rdz=1./zsl
524 cxchl=excml*rdz
525 cxchs=excms*rdz
526!
527 btgx=g/thlow
528 elfc=vkarman*btgx
529!
530 IF(pblh>1000.)THEN
531 btgh=btgx*pblh
532 ELSE
533 btgh=btgx*1000.
534 ENDIF
535!
536 wght=0.
537 wghtt=0.
538 wghtq=0.
539!----------------------------------------------------------------------
540!
541!*** SEA POINTS
542!
543!----------------------------------------------------------------------
544!
545 IF(seamask>0.5)THEN
546!
547!----------------------------------------------------------------------
548 DO itr=1,itrmx
549!----------------------------------------------------------------------
550 z0=max(ustfc*ustar*ustar,1.59e-5)
551!
552!*** VISCOUS SUBLAYER, JANJIC MWR 1994
553!
554!----------------------------------------------------------------------
555 IF(ustar<ustc)THEN
556!----------------------------------------------------------------------
557!
558 IF(ustar<ustr)THEN
559!
560 IF(ntsd==0) then
561!tgs IF(NTSD==1)THEN
562 akms=cxchs
563 akhs=cxchs
564 qs=qlow
565 ENDIF
566!
567 zu=fzu1*sqrt(sqrt(z0*ustar*rvisc))/ustar
568 wght=akms*zu*rvisc
569 rwgh=wght/(wght+1.)
570 uz0=(ulow*rwgh+uz0)*0.5
571 vz0=(vlow*rwgh+vz0)*0.5
572!
573 zt=fzt1*zu
574 zq=fzq1*zt
575 wghtt=akhs*zt*rtvisc
576 wghtq=akhs*zq*rqvisc
577!
578 IF(ntsd>0)THEN
579 thz0=((wghtt*thlow+ths)/(wghtt+1.)+thz0)*0.5
580 qz0=((wghtq*qlow+qs)/(wghtq+1.)+qz0)*0.5
581 ELSE
582 thz0=(wghtt*thlow+ths)/(wghtt+1.)
583 qz0=(wghtq*qlow+qs)/(wghtq+1.)
584 ENDIF
585!
586 ENDIF
587!
588 IF(ustar>=ustr.AND.ustar<ustc)THEN
589 zu=z0
590 uz0=0.
591 vz0=0.
592!
593 zt=fzt2*sqrt(sqrt(z0*ustar*rvisc))/ustar
594 zq=fzq2*zt
595 wghtt=akhs*zt*rtvisc
596 wghtq=akhs*zq*rqvisc
597!
598 IF(ntsd>0)THEN
599 thz0=((wghtt*thlow+ths)/(wghtt+1.)+thz0)*0.5
600 qz0=((wghtq*qlow+qs)/(wghtq+1.)+qz0)*0.5
601 ELSE
602 thz0=(wghtt*thlow+ths)/(wghtt+1.)
603 qz0=(wghtq*qlow+qs)/(wghtq+1.)
604 ENDIF
605!
606 ENDIF
607!----------------------------------------------------------------------
608 ELSE
609!----------------------------------------------------------------------
610 zu=z0
611 uz0=0.
612 vz0=0.
613!
614 zt=z0
615 thz0=ths
616!
617 zq=z0
618 qz0=qs
619!----------------------------------------------------------------------
620 ENDIF
621!----------------------------------------------------------------------
622 tem=(tlow+tz0)*0.5
623 thm=(thelow+thz0)*0.5
624!
625 a=thm*p608
626 b=(elocp/tem-1.-p608)*thm
627!
628 dthv=((thelow-thz0)*((qlow+qz0+cwmlow)*(0.5*p608)+1.) &
629 & +(qlow-qz0+cwmlow)*a+cwmlow*b)
630!
631 du2=max((ulow-uz0)**2+(vlow-vz0)**2,epsu2)
632 rib=btgx*dthv*zsl/du2
633!----------------------------------------------------------------------
634! IF(RIB>=RIC)THEN
635!----------------------------------------------------------------------
636! AKMS=MAX( VISC*RDZ,CXCHS)
637! AKHS=MAX(TVISC*RDZ,CXCHS)
638!----------------------------------------------------------------------
639! ELSE ! turbulent branch
640!----------------------------------------------------------------------
641 zslu=zsl+zu
642 zslt=zsl+zt
643!
644 rzsu=zslu/zu
645 rzst=zslt/zt
646!
647 rlogu=log(rzsu)
648 rlogt=log(rzst)
649!
650!----------------------------------------------------------------------
651!*** 1./MONIN-OBUKHOV LENGTH
652!----------------------------------------------------------------------
653!
654 rlmo=elfc*akhs*dthv/ustar**3
655!
656 zetalu=zslu*rlmo
657 zetalt=zslt*rlmo
658 zetau=zu*rlmo
659 zetat=zt*rlmo
660!
661 zetalu=min(max(zetalu,ztmin1),ztmax1)
662 zetalt=min(max(zetalt,ztmin1),ztmax1)
663 zetau=min(max(zetau,ztmin1/rzsu),ztmax1/rzsu)
664 zetat=min(max(zetat,ztmin1/rzst),ztmax1/rzst)
665!
666!----------------------------------------------------------------------
667!*** WATER FUNCTIONS
668!----------------------------------------------------------------------
669!
670 rz=(zetau-ztmin1)/dzeta1
671 k=int(rz)
672 rdzt=rz-real(k)
673 k=min(k,kztm2)
674 k=max(k,0)
675 psmz=(psim1(k+2)-psim1(k+1))*rdzt+psim1(k+1)
676!
677 rz=(zetalu-ztmin1)/dzeta1
678 k=int(rz)
679 rdzt=rz-real(k)
680 k=min(k,kztm2)
681 k=max(k,0)
682 psmzl=(psim1(k+2)-psim1(k+1))*rdzt+psim1(k+1)
683!
684 simm=psmzl-psmz+rlogu
685!
686 rz=(zetat-ztmin1)/dzeta1
687 k=int(rz)
688 rdzt=rz-real(k)
689 k=min(k,kztm2)
690 k=max(k,0)
691 pshz=(psih1(k+2)-psih1(k+1))*rdzt+psih1(k+1)
692!
693 rz=(zetalt-ztmin1)/dzeta1
694 k=int(rz)
695 rdzt=rz-real(k)
696 k=min(k,kztm2)
697 k=max(k,0)
698 pshzl=(psih1(k+2)-psih1(k+1))*rdzt+psih1(k+1)
699!
700 simh=(pshzl-pshz+rlogt)*fh01
701!----------------------------------------------------------------------
702 ustark=ustar*vkarman
703 if(abs(simm)<1.e-10.or.abs(simh)<1.e-10)then
704 write(0,*)' simm=',simm,' simh=',simh,' at i=',i,' j=',j
705 endif
706
707! if(abs(SIMM).lt.1.e-5.or.abs(SIMM).gt.1.e5)then
708 if(abs(simm).lt.1.e-10.or.abs(simh).lt.1.e-10)then
709 print*,'SIMM,PSMZL,PSMZ,RLOGU=',simm,psmzl,psmz,rlogu
710 print*,'SIMH,PSHZL,PSHZ,RLOGT,FH01=',simh,pshzl,pshz,rlogt,fh01
711 print*,'USTARK,CXCHS=',ustark,cxchs
712 print*,'PSIM1(1,2),K=',psim1(k+1),psim1(k+2),k
713 print*,'ZETAU,ZTMIN1,DZETA1=',zetau,ztmin1,dzeta1
714 print*,'PSIH1(1,2),RDZT=',psih1(k+1),psih1(k+2),rdzt
715 print*,'ZSLU,ZSLT,RLMO,ZU,ZT=',zslu,zslt,rlmo,zu,zt
716 print*,'A,B,DTHV,DU2,RIB=',a,b,dthv,du2,rib
717 errflg = 1
718 errmsg = 'ERROR(SFCDIF): '
719 return
720 end if
721
722
723
724 akms=max(ustark/simm,cxchs)
725 akhs=max(ustark/simh,cxchs)
726!
727!----------------------------------------------------------------------
728!*** BELJAARS CORRECTION FOR USTAR
729!----------------------------------------------------------------------
730!
731 IF(dthv<=0.)THEN !zj
732 wstar2=wwst2*abs(btgh*akhs*dthv)**(2./3.) !zj
733 ELSE !zj
734 wstar2=0. !zj
735 ENDIF !zj
736 ustar=max(sqrt(akms*sqrt(du2+wstar2)),epsust)
737!
738!----------------------------------------------------------------------
739! ENDIF ! End of turbulent branch
740!----------------------------------------------------------------------
741!
742 ENDDO ! End of the iteration loop over sea points
743!
744!----------------------------------------------------------------------
745!
746!*** LAND POINTS
747!
748!----------------------------------------------------------------------
749!
750 ELSE
751!
752!----------------------------------------------------------------------
753 IF(ntsd==0)THEN
754 qs=qlow
755 ENDIF
756!
757 zu=z0
758 uz0=0.
759 vz0=0.
760!
761 zt=zu*ztfc
762 thz0=ths
763!
764 zq=zt
765 qz0=qs
766!----------------------------------------------------------------------
767 tem=(tlow+tz0)*0.5
768 thm=(thelow+thz0)*0.5
769!
770 a=thm*p608
771 b=(elocp/tem-1.-p608)*thm
772!
773 dthv=((thelow-thz0)*((qlow+qz0+cwmlow)*(0.5*p608)+1.) &
774 & +(qlow-qz0+cwmlow)*a+cwmlow*b)
775!
776 du2=max((ulow-uz0)**2+(vlow-vz0)**2,epsu2)
777 rib=btgx*dthv*zsl/du2
778!----------------------------------------------------------------------
779! IF(RIB>=RIC)THEN
780! AKMS=MAX( VISC*RDZ,CXCHL)
781! AKHS=MAX(TVISC*RDZ,CXCHL)
782!----------------------------------------------------------------------
783! ELSE ! Turbulent branch
784!----------------------------------------------------------------------
785 zslu=zsl+zu
786!
787 rzsu=zslu/zu
788!
789 rlogu=log(rzsu)
790
791 zslt=zsl+zu ! u,v and t are at the same level
792!----------------------------------------------------------------------
793!
794!
795!mp Remove Topo modification of ZILFC term
796!
797! TOPOTERM=TOPOFAC*ZSFC**2.
798! TOPOTERM=MAX(TOPOTERM,3.0)
799!
800!vcw
801! RIB modification to ZILFC term
802! 7/29/2009 V Wong recommends 5, change pending
803!
804 czetmax = 10.
805! stable
806 IF(dthv>0.)THEN
807 IF (rib<ric) THEN
808 zzil=zilfc*(1.0+(rib/ric)*(rib/ric)*czetmax)
809 ELSE
810 zzil=zilfc*(czetmax+1.0)
811 ENDIF
812! unstable
813 ELSE
814 zzil=zilfc
815 ENDIF
816
817!----------------------------------------------------------------------
818!
819 land_point_iteration: DO itr=1,itrmx
820!
821!----------------------------------------------------------------------
822!*** ZILITINKEVITCH FIX FOR ZT
823!----------------------------------------------------------------------
824!
825 zt=max(exp(zzil*sqrt(ustar*zu))*zu,epszt)
826! new form ZT=MAX(EXP(ZZIL*SQRT(USTAR*Z0BASE))*Z0BASE,EPSZT)
827 rzst=zslt/zt
828 rlogt=log(rzst)
829!
830!----------------------------------------------------------------------
831!*** 1./MONIN-OBUKHOV LENGTH-SCALE
832!----------------------------------------------------------------------
833!
834 rlmo=elfc*akhs*dthv/ustar**3
835 zetalu=zslu*rlmo
836 zetalt=zslt*rlmo
837 zetau=zu*rlmo
838 zetat=zt*rlmo
839!
840 zetalu=min(max(zetalu,ztmin2),ztmax2)
841 zetalt=min(max(zetalt,ztmin2),ztmax2)
842 zetau=min(max(zetau,ztmin2/rzsu),ztmax2/rzsu)
843 zetat=min(max(zetat,ztmin2/rzst),ztmax2/rzst)
844!
845!----------------------------------------------------------------------
846!*** LAND FUNCTIONS
847!----------------------------------------------------------------------
848!
849 rz=(zetau-ztmin2)/dzeta2
850 k=int(rz)
851 rdzt=rz-real(k)
852 k=min(k,kztm2)
853 k=max(k,0)
854 psmz=(psim2(k+2)-psim2(k+1))*rdzt+psim2(k+1)
855!
856 rz=(zetalu-ztmin2)/dzeta2
857 k=int(rz)
858 rdzt=rz-real(k)
859 k=min(k,kztm2)
860 k=max(k,0)
861 psmzl=(psim2(k+2)-psim2(k+1))*rdzt+psim2(k+1)
862!
863 simm=psmzl-psmz+rlogu
864
865
866! if(abs(SIMM).lt.1.e-5.or.abs(SIMM).gt.1.e5)then
867! print*,'PSMZL,PSMZ,RLOGU=',PSMZL,PSMZ,RLOGU
868! print*,'RDZT,PSIM2(K+2)=',RDZT,PSIM2(K+2)
869! print*,'ZSL,ZU,Z0,ZSLU,ZSLT,ZT,RLMO=',ZSL,ZU,Z0,ZSLU,ZSLT,ZT,RLMO
870! print*,'RZ,K,ZETAU,ZTMIN2,DZETA2=',RZ,K,ZETAU,ZTMIN2,DZETA2
871! print*,'ELFC,AKHS,DTHV,USTAR=',ELFC,AKHS,DTHV,USTAR
872! stop
873! end if
874
875
876
877
878!
879 rz=(zetat-ztmin2)/dzeta2
880 k=int(rz)
881 rdzt=rz-real(k)
882 k=min(k,kztm2)
883 k=max(k,0)
884 pshz=(psih2(k+2)-psih2(k+1))*rdzt+psih2(k+1)
885!
886 rz=(zetalt-ztmin2)/dzeta2
887 k=int(rz)
888 rdzt=rz-real(k)
889 k=min(k,kztm2)
890 k=max(k,0)
891 pshzl=(psih2(k+2)-psih2(k+1))*rdzt+psih2(k+1)
892!
893 simh=(pshzl-pshz+rlogt)*fh02
894!----------------------------------------------------------------------
895 ustark=ustar*vkarman
896 akms=max(ustark/simm,cxchl)
897 akhs=max(ustark/simh,cxchl)
898!
899!----------------------------------------------------------------------
900!*** BELJAARS CORRECTION FOR USTAR
901!----------------------------------------------------------------------
902!
903 IF(dthv<=0.)THEN !zj
904 wstar2=wwst2*abs(btgh*akhs*dthv)**(2./3.) !zj
905 ELSE !zj
906 wstar2=0. !zj
907 ENDIF !zj
908!
909 ustar=max(sqrt(akms*sqrt(du2+wstar2)),epsust)
910!
911!----------------------------------------------------------------------
912 ENDDO land_point_iteration
913!----------------------------------------------------------------------
914!
915! ENDIF ! End of turbulant branch over land
916!
917!----------------------------------------------------------------------
918!
919 ENDIF ! End of land/sea branch
920!
921!----------------------------------------------------------------------
922!*** COUNTERGRADIENT FIX
923!----------------------------------------------------------------------
924!
925! HV=-AKHS*DTHV
926! IF(HV>0.)THEN
927! FCT=-10.*(BTGX)**(-1./3.)
928! CT=FCT*(HV/(PBLH*PBLH))**(2./3.)
929! ELSE
930 ct=0.
931! ENDIF
932!
933!----------------------------------------------------------------------
934 a1u=wght
935 a1t=wghtt
936 a1q=wghtq
937 ffm=simm*(1.+wght)
938 ffh=simh*(1.+wghtt)
939! FFQ=SIMH*(1.+WGHTQ)
940!----------------------------------------------------------------------
941!*** THE FOLLOWING DIAGNOSTIC BLOCK PRODUCES 2-m and 10-m VALUES
942!*** FOR TEMPERATURE, MOISTURE, AND WINDS. IT IS DONE HERE SINCE
943!*** THE VARIOUS QUANTITIES NEEDED FOR THE COMPUTATION ARE LOST
944!*** UPON EXIT FROM THE ROTUINE.
945!----------------------------------------------------------------------
946!----------------------------------------------------------------------
947!
948 wstar=sqrt(wstar2)/wwst
949!
950 umflx=akms*(ulow -uz0 )
951 vmflx=akms*(vlow -vz0 )
952 hsflx=akhs*(thlow-thz0)
953 hlflx=akhs*(qlow -qz0 )
954!----------------------------------------------------------------------
955! IF(RIB>=RIC)THEN
956!----------------------------------------------------------------------
957! IF(SEAMASK>0.5)THEN
958! AKMS10=MAX( VISC/10.,CXCHS)
959! AKHS02=MAX(TVISC/02.,CXCHS)
960! AKHS10=MAX(TVISC/10.,CXCHS)
961! ELSE
962! AKMS10=MAX( VISC/10.,CXCHL)
963! AKHS02=MAX(TVISC/02.,CXCHL)
964! AKHS10=MAX(TVISC/10.,CXCHL)
965! ENDIF
966!----------------------------------------------------------------------
967! ELSE
968!----------------------------------------------------------------------
969 zu10=zu+10.
970 zt02=zt+02.
971 zt10=zt+10.
972!
973 rlnu10=log(zu10/zu)
974 rlnt02=log(zt02/zt)
975 rlnt10=log(zt10/zt)
976!
977 ztau10=zu10*rlmo
978 ztat02=zt02*rlmo
979 ztat10=zt10*rlmo
980!
981!----------------------------------------------------------------------
982!*** SEA
983!----------------------------------------------------------------------
984!
985 IF(seamask>0.5)THEN
986!
987!----------------------------------------------------------------------
988 ztau10=min(max(ztau10,ztmin1),ztmax1)
989 ztat02=min(max(ztat02,ztmin1),ztmax1)
990 ztat10=min(max(ztat10,ztmin1),ztmax1)
991!----------------------------------------------------------------------
992 rz=(ztau10-ztmin1)/dzeta1
993 k=int(rz)
994 rdzt=rz-real(k)
995 k=min(k,kztm2)
996 k=max(k,0)
997 psm10=(psim1(k+2)-psim1(k+1))*rdzt+psim1(k+1)
998!
999 simm10=psm10-psmz+rlnu10
1000!
1001 rz=(ztat02-ztmin1)/dzeta1
1002 k=int(rz)
1003 rdzt=rz-real(k)
1004 k=min(k,kztm2)
1005 k=max(k,0)
1006 psh02=(psih1(k+2)-psih1(k+1))*rdzt+psih1(k+1)
1007!
1008 simh02=(psh02-pshz+rlnt02)*fh01
1009!
1010 rz=(ztat10-ztmin1)/dzeta1
1011 k=int(rz)
1012 rdzt=rz-real(k)
1013 k=min(k,kztm2)
1014 k=max(k,0)
1015 psh10=(psih1(k+2)-psih1(k+1))*rdzt+psih1(k+1)
1016!
1017 simh10=(psh10-pshz+rlnt10)*fh01
1018!
1019 akms10=max(ustark/simm10,cxchs)
1020 akhs02=max(ustark/simh02,cxchs)
1021 akhs10=max(ustark/simh10,cxchs)
1022!
1023!----------------------------------------------------------------------
1024!*** LAND
1025!----------------------------------------------------------------------
1026!
1027 ELSE
1028!
1029!----------------------------------------------------------------------
1030 ztau10=min(max(ztau10,ztmin2),ztmax2)
1031 ztat02=min(max(ztat02,ztmin2),ztmax2)
1032 ztat10=min(max(ztat10,ztmin2),ztmax2)
1033!----------------------------------------------------------------------
1034 rz=(ztau10-ztmin2)/dzeta2
1035 k=int(rz)
1036 rdzt=rz-real(k)
1037 k=min(k,kztm2)
1038 k=max(k,0)
1039 psm10=(psim2(k+2)-psim2(k+1))*rdzt+psim2(k+1)
1040!
1041 simm10=psm10-psmz+rlnu10
1042!
1043 rz=(ztat02-ztmin2)/dzeta2
1044 k=int(rz)
1045 rdzt=rz-real(k)
1046 k=min(k,kztm2)
1047 k=max(k,0)
1048 psh02=(psih2(k+2)-psih2(k+1))*rdzt+psih2(k+1)
1049!
1050 simh02=(psh02-pshz+rlnt02)*fh02
1051!
1052 rz=(ztat10-ztmin2)/dzeta2
1053 k=int(rz)
1054 rdzt=rz-real(k)
1055 k=min(k,kztm2)
1056 k=max(k,0)
1057 psh10=(psih2(k+2)-psih2(k+1))*rdzt+psih2(k+1)
1058!
1059 simh10=(psh10-pshz+rlnt10)*fh02
1060!
1061 akms10=ustark/simm10
1062 akhs02=ustark/simh02
1063 akhs10=ustark/simh10
1064!
1065 IF(akms10<=cxchl) akms10=akms
1066 IF(akhs02<=cxchl) akhs02=akhs
1067 IF(akhs10<=cxchl) akhs10=akhs
1068!
1069!----------------------------------------------------------------------
1070 ENDIF
1071!----------------------------------------------------------------------
1072! ENDIF
1073!----------------------------------------------------------------------
1074!
1075 u10 =umflx/akms10+uz0
1076 v10 =vmflx/akms10+vz0
1077 th02=hsflx/akhs02+thz0
1078 th10=hsflx/akhs10+thz0
1079 q02 =hlflx/akhs02+qz0
1080 q10 =hlflx/akhs10+qz0
1081 term1=-0.068283/tlow
1082 pshltr=psfc*exp(term1)
1083!
1084!----------------------------------------------------------------------
1085!*** COMPUTE "EQUIVALENT" Z0 TO APPROXIMATE LOCAL SHELTER READINGS.
1086!----------------------------------------------------------------------
1087!
1088 u10e=u10
1089 v10e=v10
1090!
1091 IF(seamask<0.5)THEN
1092
1093!1st ZUUZ=MIN(0.5*ZU,0.1)
1094!1st ZU=MAX(0.1*ZU,ZUUZ)
1095!tst ZUUZ=amin1(ZU*0.50,0.3)
1096!tst ZU=amax1(ZU*0.3,ZUUZ)
1097
1098 zuuz=amin1(zu*0.50,0.18)
1099 zu=amax1(zu*0.35,zuuz)
1100!
1101 zu10=zu+10.
1102 rzsu=zu10/zu
1103 rlnu10=log(rzsu)
1104
1105 zetau=zu*rlmo
1106 ztau10=zu10*rlmo
1107
1108 ztau10=min(max(ztau10,ztmin2),ztmax2)
1109 zetau=min(max(zetau,ztmin2/rzsu),ztmax2/rzsu)
1110
1111 rz=(ztau10-ztmin2)/dzeta2
1112 k=int(rz)
1113 rdzt=rz-real(k)
1114 k=min(k,kztm2)
1115 k=max(k,0)
1116 psm10=(psim2(k+2)-psim2(k+1))*rdzt+psim2(k+1)
1117 simm10=psm10-psmz+rlnu10
1118 ekms10=max(ustark/simm10,cxchl)
1119
1120 u10e=umflx/ekms10+uz0
1121 v10e=vmflx/ekms10+vz0
1122
1123 ENDIF
1124!
1125 u10=u10e
1126 v10=v10e
1127!
1128!----------------------------------------------------------------------
1129!*** SET OTHER WRF DRIVER ARRAYS
1130!----------------------------------------------------------------------
1131!
1132 rlow=plow/(r_d*tlow)
1133 chs=akhs
1134 chs2=akhs02
1135 cqs2=akhs02
1136 hfx=-rlow*cp*hsflx
1137 qfx=-rlow*hlflx*wetm
1138 flx_lh=xlv*qfx
1139 flhc=rlow*cp*akhs
1140 flqc=rlow*akhs*wetm
1141!!! QGH=PQ0/PSHLTR*EXP(A2S*(TSK-A3S)/(TSK-A4S))
1142 qgh=((1.-seamask)*pq0+seamask*pq0sea) &
1143 & /plow*exp(a2s*(tlow-a3s)/(tlow-a4s))
1144 cpm=cp*(1.+0.8*qlow)
1145!
1146!*** DO NOT COMPUTE QS OVER LAND POINTS HERE SINCE IT IS
1147!*** A PROGNOSTIC VARIABLE THERE. IT IS OKAY TO USE IT
1148!*** AS A DIAGNOSTIC OVER WATER SINCE IT WILL CAUSE NO
1149!*** INTERFERENCE BEFORE BEING RECOMPUTED IN MYJPBL.
1150!
1151 IF(seamask>0.5)THEN
1152 qs=qlow+qfx/(rlow*akhs)
1153 qs=qs/(1.-qs)
1154 ENDIF
1155!----------------------------------------------------------------------
1156 fm10=simm10+wght*simm
1157 fh2=simh02+wghtt*simh
1158!
1159 END SUBROUTINE sfcdif
1160!
1161!----------------------------------------------------------------------
1162 SUBROUTINE jsfc_init(USTAR &
1163 & ,RESTART &
1164 & ,IDS,IDE,JDS,JDE,KDS,KDE &
1165 & ,IMS,IME,JMS,JME,KMS,KME &
1166 & ,ITS,ITE,JTS,JTE,KTS,LM)
1167!----------------------------------------------------------------------
1168 IMPLICIT NONE
1169!----------------------------------------------------------------------
1170 LOGICAL,INTENT(IN) :: restart
1171!
1172 INTEGER,INTENT(IN) :: ids,ide,jds,jde,kds,kde &
1173 & ,IMS,IME,JMS,JME,KMS,KME &
1174 & ,ITS,ITE,JTS,JTE,KTS,LM
1175!
1176 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: ustar
1177!
1178!----------------------------------------------------------------------
1179!*** LOCAL VARIABLES
1180!----------------------------------------------------------------------
1181 INTEGER :: i,j,k,itf,jtf,ktf
1182!
1183 REAL(kind=kfpt) :: x,zeta1,zeta2,zrng1,zrng2
1184!
1185 REAL(kind=kfpt) :: pihf=3.1415926/2.,eps=1.e-6
1186!
1187!----------------------------------------------------------------------
1188 jtf=min0(jte,jde-1)
1189 ktf=min0(lm,kde-1)
1190 itf=min0(ite,ide-1)
1191!
1192!
1193!*** FOR NOW, ASSUME SIGMA MODE FOR LOWEST MODEL LAYER
1194!
1195!----------------------------------------------------------------------
1196 IF(.NOT.restart)THEN
1197 DO j=jts,jte
1198 DO i=its,itf
1199 ustar(i,j)=0.1
1200 ENDDO
1201 ENDDO
1202 ENDIF
1203
1204!----------------------------------------------------------------------
1205!
1206!*** COMPUTE SURFACE LAYER INTEGRAL FUNCTIONS
1207!
1208!----------------------------------------------------------------------
1209 fh01=1.
1210 fh02=1.
1211!
1212! ZTMIN1=-10.0
1213! ZTMAX1=2.0
1214! ZTMIN2=-10.0
1215! ZTMAX2=2.0
1216!org b
1217! ZTMIN1=-5.0
1218! ZTMAX1=1.0
1219! ZTMIN2=-5.0
1220! ZTMAX2=1.0
1221!org e
1222 ztmin1=-5.0
1223 ztmax1=9.0
1224 ztmin2=-5.0
1225 ztmax2=9.0
1226
1227 zrng1=ztmax1-ztmin1
1228 zrng2=ztmax2-ztmin2
1229!
1230 dzeta1=zrng1/(kztm-1)
1231 dzeta2=zrng2/(kztm-1)
1232!
1233!----------------------------------------------------------------------
1234!*** FUNCTION DEFINITION LOOP
1235!----------------------------------------------------------------------
1236!
1237 zeta1=ztmin1
1238 zeta2=ztmin2
1239!
1240 DO k=1,kztm
1241!
1242!----------------------------------------------------------------------
1243!*** UNSTABLE RANGE
1244!----------------------------------------------------------------------
1245!
1246 IF(zeta1<0.)THEN
1247!
1248!----------------------------------------------------------------------
1249!*** PAULSON 1970 FUNCTIONS
1250!----------------------------------------------------------------------
1251 x=sqrt(sqrt(1.-16.*zeta1))
1252!
1253 psim1(k)=-2.*log((x+1.)/2.)-log((x*x+1.)/2.)+2.*atan(x)-pihf
1254 psih1(k)=-2.*log((x*x+1.)/2.)
1255!
1256!----------------------------------------------------------------------
1257!*** STABLE RANGE
1258!----------------------------------------------------------------------
1259!
1260 ELSE
1261!
1262!----------------------------------------------------------------------
1263!*** PAULSON 1970 FUNCTIONS
1264!----------------------------------------------------------------------
1265!
1266! PSIM1(K)=5.*ZETA1
1267! PSIH1(K)=5.*ZETA1
1268!----------------------------------------------------------------------
1269!*** HOLTSLAG AND DE BRUIN 1988
1270!----------------------------------------------------------------------
1271!
1272 psim1(k)=0.7*zeta1+0.75*zeta1*(6.-0.35*zeta1)*exp(-0.35*zeta1)
1273 psih1(k)=0.7*zeta1+0.75*zeta1*(6.-0.35*zeta1)*exp(-0.35*zeta1)
1274!----------------------------------------------------------------------
1275!
1276 ENDIF
1277!
1278!----------------------------------------------------------------------
1279!*** UNSTABLE RANGE
1280!----------------------------------------------------------------------
1281!
1282 IF(zeta2<0.)THEN
1283!
1284!----------------------------------------------------------------------
1285!*** PAULSON 1970 FUNCTIONS
1286!----------------------------------------------------------------------
1287!
1288 x=sqrt(sqrt(1.-16.*zeta2))
1289!
1290 psim2(k)=-2.*log((x+1.)/2.)-log((x*x+1.)/2.)+2.*atan(x)-pihf
1291 psih2(k)=-2.*log((x*x+1.)/2.)
1292!----------------------------------------------------------------------
1293!*** STABLE RANGE
1294!----------------------------------------------------------------------
1295!
1296 ELSE
1297!
1298!----------------------------------------------------------------------
1299!*** PAULSON 1970 FUNCTIONS
1300!----------------------------------------------------------------------
1301!
1302! PSIM2(K)=5.*ZETA2
1303! PSIH2(K)=5.*ZETA2
1304!
1305!----------------------------------------------------------------------
1306!*** HOLTSLAG AND DE BRUIN 1988
1307!----------------------------------------------------------------------
1308!
1309 psim2(k)=0.7*zeta2+0.75*zeta2*(6.-0.35*zeta2)*exp(-0.35*zeta2)
1310 psih2(k)=0.7*zeta2+0.75*zeta2*(6.-0.35*zeta2)*exp(-0.35*zeta2)
1311!----------------------------------------------------------------------
1312!
1313 ENDIF
1314!
1315!----------------------------------------------------------------------
1316 IF(k==kztm)THEN
1317 ztmax1=zeta1
1318 ztmax2=zeta2
1319 ENDIF
1320!
1321 zeta1=zeta1+dzeta1
1322 zeta2=zeta2+dzeta2
1323!----------------------------------------------------------------------
1324 ENDDO
1325!----------------------------------------------------------------------
1326 ztmax1=ztmax1-eps
1327 ztmax2=ztmax2-eps
1328!----------------------------------------------------------------------
1329!
1330 END SUBROUTINE jsfc_init
1331!
1332!----------------------------------------------------------------------
1333!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
1334!-----------------------------------------------------------------------
1335!
1336 END MODULE module_sf_jsfc
1337!
1338!-----------------------------------------------------------------------
subroutine sfcdif
This subroutine calculates surface layer exchange coefficients via iterative process(see Chen et al....
Definition sflx.f:1989