source: lmdz_wrf/WRFV3/phys/module_sf_myjsfc.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

  • Property svn:executable set to *
File size: 42.8 KB
Line 
1!----------------------------------------------------------------------
2!
3      MODULE MODULE_SF_MYJSFC
4!
5!----------------------------------------------------------------------
6!
7      USE MODULE_MODEL_CONSTANTS
8#ifdef DM_PARALLEL
9      USE MODULE_DM, ONLY : WRF_DM_MAXVAL
10#endif
11!
12!----------------------------------------------------------------------
13!
14! REFERENCES:  Janjic (2002), NCEP Office Note 437
15!
16! ABSTRACT:
17!     MYJSFC GENERATES THE SURFACE EXCHANGE COEFFICIENTS FOR VERTICAL
18!     TURBULENT EXCHANGE BASED UPON MONIN_OBUKHOV THEORY WITH
19!     VARIOUS REFINEMENTS.
20!
21!----------------------------------------------------------------------
22!
23      INTEGER :: ITRMX=5 ! Iteration count for sfc layer computations
24!
25      REAL,PARAMETER :: VKARMAN=0.4
26      REAL,PARAMETER :: CAPA=R_D/CP,ELOCP=2.72E6/CP,RCAP=1./CAPA
27      REAL,PARAMETER :: GOCP02=G/CP*2.,GOCP10=G/CP*10.
28!      REAL,PARAMETER :: EPSU2=1.E-4,EPSUST=0.07,EPSZT=1.E-5
29!       ECMWF sets lower limit on ln(Z0h)=-20 --> EPSZT=2.e-9
30      REAL,PARAMETER :: EPSU2=1.E-6,EPSUST=1.e-9,EPSZT=1.E-28
31      REAL,PARAMETER :: A2S=17.2693882,A3S=273.16,A4S=35.86
32      REAL,PARAMETER :: SEAFC=0.98,PQ0SEA=PQ0*SEAFC
33      REAL,PARAMETER :: BETA=1./273.,EXCML=0.0001,EXCMS=0.0001         &
34     &                 ,GLKBR=10.,GLKBS=30.,PI=3.1415926               &
35     &                 ,QVISC=2.1E-5,RIC=0.505,SMALL=0.35              &
36     &                 ,SQPR=0.84,SQSC=0.84,SQVISC=258.2,TVISC=2.1E-5  &
37     &                 ,USTC=0.7,USTR=0.225,VISC=1.5E-5,FH=1.01        &
38     &                 ,WWST=1.2,ZTFC=1.,TOPOFAC=9.0e-6
39!
40      REAL,PARAMETER :: BTG=BETA*G,CZIV=SMALL*GLKBS                    &
41     &                 ,GRRS=GLKBR/GLKBS                               &
42     &                 ,RTVISC=1./TVISC,RVISC=1./VISC                  &
43     &                 ,ZQRZT=SQSC/SQPR                                &
44     &                 ,FZQ1=RTVISC*QVISC*ZQRZT                        &
45     &                 ,FZQ2=RTVISC*QVISC*ZQRZT                        &
46     &                 ,FZT1=RVISC *TVISC*SQPR                         &
47     &                 ,FZT2=CZIV*GRRS*TVISC*SQPR                      &
48     &                 ,FZU1=CZIV*VISC                                 &
49     &                 ,PIHF=0.5*PI                                    &
50     &                 ,RQVISC=1./QVISC                                &
51     &                 ,USTFC=0.018/G                                  &
52     &                 ,WWST2=WWST*WWST
53!
54!----------------------------------------------------------------------
55      INTEGER, PARAMETER :: KZTM=10001,KZTM2=KZTM-2
56!
57      REAL :: DZETA1,DZETA2,FH01,FH02,ZTMAX1,ZTMAX2,ZTMIN1,ZTMIN2
58!
59      REAL,DIMENSION(KZTM) :: PSIH1,PSIH2,PSIM1,PSIM2
60!
61!----------------------------------------------------------------------
62!
63CONTAINS
64!
65!----------------------------------------------------------------------
66      SUBROUTINE MYJSFC(ITIMESTEP,HT,DZ                                &
67     &                 ,PMID,PINT,TH,T,QV,QC,U,V,Q2                    &
68     &                 ,TSK,QSFC,THZ0,QZ0,UZ0,VZ0                      &
69     &                 ,LOWLYR,XLAND,IVGTYP,ISURBAN,IZ0TLND            &
70     &                 ,USTAR,ZNT,Z0BASE,PBLH,MAVAIL,RMOL              &
71     &                 ,AKHS,AKMS                                      &
72     &                 ,RIB                                            &
73     &                 ,CHS,CHS2,CQS2,HFX,QFX,FLX_LH,FLHC,FLQC         &
74     &                 ,QGH,CPM,CT                                     &
75     &                 ,U10,V10,T02,TH02,TSHLTR,TH10,Q02,QSHLTR,Q10,PSHLTR          &
76     &                 ,P1000mb                                        &
77     &                 ,IDS,IDE,JDS,JDE,KDS,KDE                        &
78     &                 ,IMS,IME,JMS,JME,KMS,KME                        &
79     &                 ,ITS,ITE,JTS,JTE,KTS,KTE)
80!----------------------------------------------------------------------
81!
82      IMPLICIT NONE
83!
84!----------------------------------------------------------------------
85      INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                    &
86     &                     ,IMS,IME,JMS,JME,KMS,KME                    &
87     &                     ,ITS,ITE,JTS,JTE,KTS,KTE
88!
89      INTEGER,INTENT(IN) :: ITIMESTEP
90!
91      INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: LOWLYR
92!
93      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: HT,MAVAIL,TSK      &
94     &                                             ,XLAND,Z0BASE
95!
96      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: DZ         &
97     &                                                     ,PMID,PINT  &
98     &                                                     ,Q2,QC,QV   &
99     &                                                     ,T,TH       &
100     &                                                     ,U,V
101      INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IVGTYP
102      INTEGER,                           INTENT(IN) :: ISURBAN
103      INTEGER,                           INTENT(IN) :: IZ0TLND
104!
105      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: FLX_LH,HFX,PSHLTR &
106     &                                              ,QFX,Q10,QSHLTR    &
107     &                                              ,TH10,TSHLTR,T02   &
108     &                                              ,U10,V10,TH02,Q02
109!
110      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: AKHS,AKMS       &
111     &                                                ,PBLH,QSFC
112      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT)   :: RIB
113!
114      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: QZ0,RMOL,THZ0   &
115     &                                                ,USTAR,UZ0,VZ0   &
116     &                                                ,ZNT
117!
118      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: CHS,CHS2,CQS2     &
119     &                                              ,CPM,CT,FLHC,FLQC  &
120     &                                              ,QGH
121!
122      REAL,INTENT(IN) :: P1000mb
123!----------------------------------------------------------------------
124!***
125!***  LOCAL VARIABLES
126!***
127      INTEGER :: I,J,K,KFLIP,LMH,LPBL,NTSD
128!
129      REAL :: APESFC,BTGX,CWMLOW                                       &
130     &       ,DQDT,DTDIF,DTDT,DUDT,DVDT                                &
131     &       ,FIS                                                      &
132     &       ,P02P,P10P,PLOW,PSFC,PTOP,QLOW,QS02,QS10                  &
133     &       ,RAPA,RAPA02,RAPA10,RATIOMX,RDZ,SEAMASK,SM                &
134     &       ,T02P,T10P,TH02P,TH10P,THLOW,THvLOW                       &
135     &       ,TLOW,TZ0,ULOW,VLOW,ZSL
136!
137      REAL,DIMENSION(KTS:KTE) :: CWMK,PK,Q2K,QK,THvK,THK,TK,UK,VK
138!
139      REAL,DIMENSION(KTS:KTE-1) :: EL,ELM
140!
141      REAL,DIMENSION(KTS:KTE+1) :: ZHK
142!
143      REAL,DIMENSION(ITS:ITE,JTS:JTE) :: THSK
144!
145      REAL,DIMENSION(ITS:ITE,KTS:KTE+1,JTS:JTE) :: ZINT
146!
147!----------------------------------------------------------------------
148!**********************************************************************
149!----------------------------------------------------------------------
150!
151!***  MAKE PREPARATIONS
152!
153!----------------------------------------------------------------------
154      DO J=JTS,JTE
155      DO K=KTS,KTE+1
156      DO I=ITS,ITE
157        ZINT(I,K,J)=0.
158      ENDDO
159      ENDDO
160      ENDDO
161!
162      DO J=JTS,JTE
163      DO I=ITS,ITE
164        ZINT(I,KTE+1,J)=HT(I,J)     ! Z at bottom of lowest sigma layer
165        PBLH(I,J)=-1.
166!
167!!!!!!!!!
168!!!!!! UNCOMMENT THESE LINES IF USING ETA COORDINATES
169!!!!!!!!!
170!!!!!!  ZINT(I,KTE+1,J)=1.E-4       ! Z of bottom of lowest eta layer
171!!!!!!  ZHK(KTE+1)=1.E-4            ! Z of bottom of lowest eta layer
172!
173      ENDDO
174      ENDDO
175!
176      DO J=JTS,JTE
177      DO K=KTE,KTS,-1
178        KFLIP=KTE+1-K
179        DO I=ITS,ITE
180          ZINT(I,K,J)=ZINT(I,K+1,J)+DZ(I,KFLIP,J)
181        ENDDO
182      ENDDO
183      ENDDO
184!
185      NTSD=ITIMESTEP
186!
187#if ( NMM_CORE == 1 )
188     if(NTSD+1.eq.1) then
189#else
190     IF(NTSD==1)THEN
191#endif
192!tgs       IF(NTSD==1)THEN
193        DO J=JTS,JTE
194        DO I=ITS,ITE
195          USTAR(I,J)=0.1
196          FIS=HT(I,J)*G
197          SM=XLAND(I,J)-1.
198!!!       Z0(I,J)=SM*Z0SEA+(1.-SM)*(Z0(I,J)*Z0MAX+FIS*FCM+Z0LAND)
199!!!       ZNT(I,J)=SM*Z0SEA+(1.-SM)*(ZNT(I,J)*Z0MAX+FIS*FCM+Z0LAND)
200        ENDDO
201        ENDDO
202      ENDIF
203!
204!!!!  IF(NTSD==1)THEN
205        DO J=JTS,JTE
206        DO I=ITS,ITE
207          CT(I,J)=0.
208        ENDDO
209        ENDDO
210!!!!  ENDIF
211!
212!----------------------------------------------------------------------
213        setup_integration:  DO J=JTS,JTE
214!----------------------------------------------------------------------
215!
216        DO I=ITS,ITE
217!
218!***  LOWEST LAYER ABOVE GROUND MUST BE FLIPPED
219!
220          LMH=KTE-LOWLYR(I,J)+1
221!
222          PTOP=PINT(I,KTE+1,J)      ! KTE+1=KME
223          PSFC=PINT(I,LOWLYR(I,J),J)
224! Define THSK here (for first timestep mostly)
225!          THSK(I,J)=TSK(I,J)/(PSFC*1.E-5)**CAPA
226          THSK(I,J)=TSK(I,J)/(PSFC/P1000mb)**CAPA
227!
228!***  CONVERT LAND MASK (1 FOR SEA; 0 FOR LAND)
229!
230          SEAMASK=XLAND(I,J)-1.
231!
232!***  FILL 1-D VERTICAL ARRAYS
233!***  AND FLIP DIRECTION SINCE MYJ SCHEME
234!***  COUNTS DOWNWARD FROM THE DOMAIN'S TOP
235!
236          DO K=KTE,KTS,-1
237            KFLIP=KTE+1-K
238            THK(K)=TH(I,KFLIP,J)
239            TK(K)=T(I,KFLIP,J)
240            RATIOMX=QV(I,KFLIP,J)
241            QK(K)=RATIOMX/(1.+RATIOMX)
242            PK(K)=PMID(I,KFLIP,J)
243            CWMK(K)=QC(I,KFLIP,J)
244            THvK(K)=(qk(k)*0.608+(1.-CWMK(K)))*THK(K)
245            Q2K(K)=2.*Q2(I,KFLIP,J)
246!
247!***  COMPUTE THE HEIGHTS OF THE LAYER INTERFACES
248!
249            ZHK(K)=ZINT(I,K,J)
250!
251          ENDDO
252          ZHK(KTE+1)=HT(I,J)          ! Z at bottom of lowest sigma layer
253!
254          DO K=KTE,KTS,-1
255            KFLIP=KTE+1-K
256            UK(K)=U(I,KFLIP,J)
257            VK(K)=V(I,KFLIP,J)
258          ENDDO
259!
260!***  FIND THE HEIGHT OF THE PBL
261!
262          LPBL=LMH
263          DO K=LMH-1,1,-1
264            IF(Q2K(K)<=EPSQ2*FH) THEN
265              LPBL=K
266              GO TO 110
267            ENDIF
268          ENDDO
269!
270          LPBL=1
271!
272!-----------------------------------------------------------------------
273!--------------THE HEIGHT OF THE PBL------------------------------------
274!-----------------------------------------------------------------------
275!
276 110      PBLH(I,J)=ZHK(LPBL+1)-ZHK(LMH+1)
277!
278!----------------------------------------------------------------------
279!***
280!***  FIND THE SURFACE EXCHANGE COEFFICIENTS
281!***
282!----------------------------------------------------------------------
283          PLOW=PK(LMH)
284          TLOW=TK(LMH)
285          THLOW=THK(LMH)
286          THvLOW=THvK(LMH)
287          QLOW=QK(LMH)
288          CWMLOW=CWMK(LMH)
289          ULOW=UK(LMH)
290          VLOW=VK(LMH)
291          ZSL=(ZHK(LMH)-ZHK(LMH+1))*0.5
292!          APESFC=(PSFC*1.E-5)**CAPA
293          APESFC=(PSFC/P1000mb)**CAPA
294!tgs - in ARW THZ0 is not initialized when MYJSFC is called first time
295#if ( NMM_CORE == 1 )
296     if(NTSD+1.eq.1) then
297#else
298     IF(NTSD==1)THEN
299#endif
300!       if(itimestep.le.1) then
301          TZ0=TSK(I,J)
302       else
303          TZ0=THZ0(I,J)*APESFC
304       endif
305!
306          CALL SFCDIF(NTSD,SEAMASK,THSK(I,J),QSFC(I,J),PSFC            &
307     &               ,UZ0(I,J),VZ0(I,J),TZ0,TSK(I,J),THZ0(I,J),QZ0(I,J)         &
308     &               ,USTAR(I,J),ZNT(I,J),Z0BASE(I,J),CT(I,J),RMOL(I,J) &
309     &               ,AKMS(I,J),AKHS(I,J),RIB(I,J),PBLH(I,J),MAVAIL(I,J) &
310     &               ,IVGTYP(I,J),ISURBAN,IZ0TLND                      &
311     &               ,CHS(I,J),CHS2(I,J),CQS2(I,J)                     &
312     &               ,HFX(I,J),QFX(I,J),FLX_LH(I,J)                    &
313     &               ,FLHC(I,J),FLQC(I,J),QGH(I,J),CPM(I,J)            &
314     &               ,ULOW,VLOW,TLOW,THLOW,THvLOW,QLOW,CWMLOW          &
315     &               ,ZSL,PLOW                                         &
316     &               ,U10(I,J),V10(I,J),TSHLTR(I,J),TH10(I,J)          &
317     &               ,QSHLTR(I,J),Q10(I,J),PSHLTR(I,J)                 &
318     &               ,IDS,IDE,JDS,JDE,KDS,KDE                          &
319     &               ,IMS,IME,JMS,JME,KMS,KME                          &
320     &               ,ITS,ITE,JTS,JTE,KTS,KTE,i,j,ZHK(LMH+1))
321!
322!***  REMOVE SUPERATURATION AT 2M AND 10M
323!
324          RAPA=APESFC
325          TH02P=TSHLTR(I,J)
326          TH10P=TH10(I,J)
327          TH02(I,J)=TSHLTR(I,J)
328!
329          RAPA02=RAPA-GOCP02/TH02P
330          RAPA10=RAPA-GOCP10/TH10P
331!
332          T02P=TH02P*RAPA02
333          T10P=TH10P*RAPA10
334! 1 may 06 tgs         T02(I,J) = T02P
335          T02(I,J) = TH02(I,J)*APESFC
336!
337!          P02P=(RAPA02**RCAP)*1.E5
338!          P10P=(RAPA10**RCAP)*1.E5
339          P02P=(RAPA02**RCAP)*P1000mb
340          P10P=(RAPA10**RCAP)*P1000mb
341!
342          QS02=PQ0/P02P*EXP(A2*(T02P-A3)/(T02P-A4))
343          QS10=PQ0/P10P*EXP(A2*(T10P-A3)/(T10P-A4))
344!
345          IF(QSHLTR(I,J)>QS02)QSHLTR(I,J)=QS02
346          IF(Q10   (I,J)>QS10)Q10   (I,J)=QS10
347          Q02(I,J)=QSHLTR(I,J)/(1.-QSHLTR(I,J))
348!----------------------------------------------------------------------
349!
350        ENDDO
351!
352!----------------------------------------------------------------------
353      ENDDO setup_integration
354!----------------------------------------------------------------------
355
356      END SUBROUTINE MYJSFC
357!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
358!----------------------------------------------------------------------
359      SUBROUTINE SFCDIF(NTSD,SEAMASK,THS,QS,PSFC                       &
360     &                 ,UZ0,VZ0,TZ0,TSK,THZ0,QZ0                       &
361     &                 ,USTAR,Z0,Z0BASE,CT,RLMO,AKMS,AKHS,RIB,PBLH,WETM &
362     &                 ,VEGTYP,ISURBAN,IZ0TLND                         &
363     &                 ,CHS,CHS2,CQS2,HFX,QFX,FLX_LH,FLHC,FLQC,QGH,CPM &
364     &                 ,ULOW,VLOW,TLOW,THLOW,THvLOW,QLOW,CWMLOW        &
365     &                 ,ZSL,PLOW                                       &
366     &                 ,U10,V10,TH02,TH10,Q02,Q10,PSHLTR               &
367     &                 ,IDS,IDE,JDS,JDE,KDS,KDE                        &
368     &                 ,IMS,IME,JMS,JME,KMS,KME                        &
369     &                 ,ITS,ITE,JTS,JTE,KTS,KTE,i,j,ZSFC)
370!     ****************************************************************
371!     *                                                              *
372!     *                       SURFACE LAYER                          *
373!     *                                                              *
374!     ****************************************************************
375!----------------------------------------------------------------------
376!
377      IMPLICIT NONE
378!
379!----------------------------------------------------------------------
380      INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                    &
381     &                     ,IMS,IME,JMS,JME,KMS,KME                    &
382     &                     ,ITS,ITE,JTS,JTE,KTS,KTE,i,j
383!
384      INTEGER,INTENT(IN) :: NTSD
385!
386      REAL,INTENT(IN) :: CWMLOW,PBLH,PLOW,QLOW,PSFC,SEAMASK,ZSFC       &
387     &                  ,THvLOW,THLOW,THS,TSK,TLOW,TZ0                 &
388     &                  ,ULOW,VLOW,WETM,ZSL                            &
389     &                  ,Z0BASE
390      INTEGER, INTENT(IN) :: VEGTYP, ISURBAN, IZ0TLND
391!
392      REAL,INTENT(OUT) :: CHS,CHS2,CPM,CQS2,CT,FLHC,FLQC,FLX_LH,HFX    &
393     &                   ,PSHLTR,Q02,Q10,QFX,QGH,RIB,RLMO              &
394     &                   ,TH02,TH10,U10,V10
395!
396      REAL,INTENT(INOUT) :: AKHS,AKMS,QZ0,THZ0,USTAR,UZ0,VZ0,Z0,QS
397!----------------------------------------------------------------------
398!***
399!***  LOCAL VARIABLES
400!***
401      INTEGER :: ITR,K
402!
403      REAL :: BTGH,BTGX,CXCHL,CXCHS,CZETMAX,DTHV,DU2,ELFC,FCT          &
404     &       ,HLFLX,HSFLX,HV,PSH02,PSH10,PSHZ,PSHZL,PSM10,PSMZ,PSMZL   &
405     &       ,RDZ,RDZT,RLMA,RLMN,RLMP                                  &
406     &       ,RLOGT,RLOGU,RWGH,RZ,RZST,RZSU,SIMH,SIMM                  &
407     &       ,UMFLX,USTARK,VMFLX,WGHT,WGHTT,WGHTQ,WSTAR2               &
408     &       ,X,XLT,XLT4,XLU,XLU4,XT,XT4,XU,XU4,ZETALT,ZETALU          &
409     &       ,ZETAT,ZETAU,ZQ,ZSLT,ZSLU,ZT,ZU,TOPOTERM,ZZIL
410!
411!***  DIAGNOSTICS
412!
413      REAL :: AKHS02,AKHS10,AKMS02,AKMS10,EKMS10,QSAT10,QSAT2          &
414     &       ,RLNT02,RLNT10,RLNU10,SIMH02,SIMH10,SIMM10,T02,T10        &
415     &       ,TERM1,RLOW,U10E,V10E,WSTAR,XLT02,XLT024,XLT10            &
416     &       ,XLT104,XLU10,XLU104,XU10,XU104,ZT02,ZT10,ZTAT02,ZTAT10   &
417     &       ,ZTAU,ZTAU10,ZU10,ZUUZ
418      REAL :: CZIL
419      REAL :: ZILFC
420!----------------------------------------------------------------------
421!**********************************************************************
422!----------------------------------------------------------------------
423      RDZ=1./ZSL
424      CXCHL=EXCML*RDZ
425      CXCHS=EXCMS*RDZ
426!
427      BTGX=G/THLOW
428      ELFC=VKARMAN*BTGX
429!
430      IF(PBLH>1000.)THEN
431        BTGH=BTGX*PBLH
432      ELSE
433        BTGH=BTGX*1000.
434      ENDIF
435!
436!----------------------------------------------------------------------
437!
438!***  SEA POINTS
439!
440!----------------------------------------------------------------------
441!
442      IF(SEAMASK>0.5)THEN
443!
444!----------------------------------------------------------------------
445        DO ITR=1,ITRMX
446!----------------------------------------------------------------------
447          Z0=MAX(USTFC*USTAR*USTAR,1.59E-5)
448!
449!***  VISCOUS SUBLAYER, JANJIC MWR 1994
450!
451!----------------------------------------------------------------------
452          IF(USTAR<USTC)THEN
453!----------------------------------------------------------------------
454!
455            IF(USTAR<USTR)THEN
456!
457#if ( NMM_CORE == 1 )
458     if(NTSD+1.eq.1) then
459#else
460     IF(NTSD==1)THEN
461#endif
462!tgs              IF(NTSD==1)THEN
463                AKMS=CXCHS
464                AKHS=CXCHS
465                QS=QLOW
466              ENDIF
467!
468              ZU=FZU1*SQRT(SQRT(Z0*USTAR*RVISC))/USTAR
469              WGHT=AKMS*ZU*RVISC
470              RWGH=WGHT/(WGHT+1.)
471              UZ0=(ULOW*RWGH+UZ0)*0.5
472              VZ0=(VLOW*RWGH+VZ0)*0.5
473!
474              ZT=FZT1*ZU
475              ZQ=FZQ1*ZT
476              WGHTT=AKHS*ZT*RTVISC
477              WGHTQ=AKHS*ZQ*RQVISC
478!
479#if ( NMM_CORE == 1 )
480     if(NTSD+1>1) then
481#else
482     IF(NTSD>1)THEN
483#endif
484!tgs              IF(NTSD>1)THEN
485                THZ0=((WGHTT*THLOW+THS)/(WGHTT+1.)+THZ0)*0.5
486                QZ0=((WGHTQ*QLOW+QS)/(WGHTQ+1.)+QZ0)*0.5
487              ELSE
488                THZ0=(WGHTT*THLOW+THS)/(WGHTT+1.)
489                QZ0=(WGHTQ*QLOW+QS)/(WGHTQ+1.)
490              ENDIF
491!
492            ENDIF
493!
494            IF(USTAR>=USTR.AND.USTAR<USTC)THEN
495              ZU=Z0
496              UZ0=0.
497              VZ0=0.
498!
499              ZT=FZT2*SQRT(SQRT(Z0*USTAR*RVISC))/USTAR
500              ZQ=FZQ2*ZT
501              WGHTT=AKHS*ZT*RTVISC
502              WGHTQ=AKHS*ZQ*RQVISC
503!
504#if ( NMM_CORE == 1 )
505     if(NTSD+1>1) then
506#else
507     IF(NTSD>1)THEN
508#endif
509
510!tgs              IF(NTSD>1)THEN
511                THZ0=((WGHTT*THLOW+THS)/(WGHTT+1.)+THZ0)*0.5
512                QZ0=((WGHTQ*QLOW+QS)/(WGHTQ+1.)+QZ0)*0.5
513              ELSE
514                THZ0=(WGHTT*THLOW+THS)/(WGHTT+1.)
515                QZ0=(WGHTQ*QLOW+QS)/(WGHTQ+1.)
516              ENDIF
517!
518            ENDIF
519!----------------------------------------------------------------------
520          ELSE
521!----------------------------------------------------------------------
522            ZU=Z0
523            UZ0=0.
524            VZ0=0.
525!
526            ZT=Z0
527            THZ0=THS
528!
529            ZQ=Z0
530            QZ0=QS
531!----------------------------------------------------------------------
532          ENDIF
533!----------------------------------------------------------------------
534          dthv=thvlow-thz0*(0.608*qz0+1.)
535!
536          DU2=MAX((ULOW-UZ0)**2+(VLOW-VZ0)**2,EPSU2)
537          RIB=BTGX*DTHV*ZSL/DU2
538!----------------------------------------------------------------------
539!         IF(RIB>=RIC)THEN
540!----------------------------------------------------------------------
541!           AKMS=MAX( VISC*RDZ,CXCHS)
542!           AKHS=MAX(TVISC*RDZ,CXCHS)
543!----------------------------------------------------------------------
544!         ELSE  !  turbulent branch
545!----------------------------------------------------------------------
546            ZSLU=ZSL+ZU
547            ZSLT=ZSL+ZT
548!
549            RZSU=ZSLU/ZU
550            RZST=ZSLT/ZT
551!
552            RLOGU=LOG(RZSU)
553            RLOGT=LOG(RZST)
554!
555!----------------------------------------------------------------------
556!***  1./MONIN-OBUKHOV LENGTH
557!----------------------------------------------------------------------
558!
559            RLMO=ELFC*AKHS*DTHV/USTAR**3
560!
561            ZETALU=ZSLU*RLMO
562            ZETALT=ZSLT*RLMO
563            ZETAU=ZU*RLMO
564            ZETAT=ZT*RLMO
565!
566            ZETALU=MIN(MAX(ZETALU,ZTMIN1),ZTMAX1)
567            ZETALT=MIN(MAX(ZETALT,ZTMIN1),ZTMAX1)
568            ZETAU=MIN(MAX(ZETAU,ZTMIN1/RZSU),ZTMAX1/RZSU)
569            ZETAT=MIN(MAX(ZETAT,ZTMIN1/RZST),ZTMAX1/RZST)
570!
571!----------------------------------------------------------------------
572!***   WATER FUNCTIONS
573!----------------------------------------------------------------------
574!
575            RZ=(ZETAU-ZTMIN1)/DZETA1
576            K=INT(RZ)
577            RDZT=RZ-REAL(K)
578            K=MIN(K,KZTM2)
579            K=MAX(K,0)
580            PSMZ=(PSIM1(K+2)-PSIM1(K+1))*RDZT+PSIM1(K+1)
581!
582            RZ=(ZETALU-ZTMIN1)/DZETA1
583            K=INT(RZ)
584            RDZT=RZ-REAL(K)
585            K=MIN(K,KZTM2)
586            K=MAX(K,0)
587            PSMZL=(PSIM1(K+2)-PSIM1(K+1))*RDZT+PSIM1(K+1)
588!
589            SIMM=PSMZL-PSMZ+RLOGU
590!
591            RZ=(ZETAT-ZTMIN1)/DZETA1
592            K=INT(RZ)
593            RDZT=RZ-REAL(K)
594            K=MIN(K,KZTM2)
595            K=MAX(K,0)
596            PSHZ=(PSIH1(K+2)-PSIH1(K+1))*RDZT+PSIH1(K+1)
597!
598            RZ=(ZETALT-ZTMIN1)/DZETA1
599            K=INT(RZ)
600            RDZT=RZ-REAL(K)
601            K=MIN(K,KZTM2)
602            K=MAX(K,0)
603            PSHZL=(PSIH1(K+2)-PSIH1(K+1))*RDZT+PSIH1(K+1)
604!
605            SIMH=(PSHZL-PSHZ+RLOGT)*FH01
606
607!----------------------------------------------------------------------
608            USTARK=USTAR*VKARMAN
609!!!!!!!!!!!!!!!!!!!!!
610            AKMS=MAX(USTARK/SIMM,CXCHS)
611            AKHS=MAX(USTARK/SIMH,CXCHS)
612!
613!----------------------------------------------------------------------
614!***  BELJAARS CORRECTION FOR USTAR
615!----------------------------------------------------------------------
616!
617            IF(DTHV<=0.)THEN                                           !zj
618              WSTAR2=WWST2*ABS(BTGH*AKHS*DTHV)**(2./3.)                !zj
619            ELSE                                                       !zj
620              WSTAR2=0.                                                !zj
621            ENDIF                                                      !zj
622            USTAR=MAX(SQRT(AKMS*SQRT(DU2+WSTAR2)),EPSUST)
623!
624!----------------------------------------------------------------------
625!         ENDIF  !  End of turbulent branch
626!----------------------------------------------------------------------
627!
628        ENDDO  !  End of the iteration loop over sea points
629!
630!----------------------------------------------------------------------
631!
632!***  LAND POINTS
633!
634!----------------------------------------------------------------------
635!
636      ELSE
637!
638!----------------------------------------------------------------------
639        IF(NTSD==1)THEN
640          QS=QLOW
641        ENDIF
642!
643        ZU=Z0
644        UZ0=0.
645        VZ0=0.
646!
647        ZT=ZU*ZTFC
648        THZ0=THS
649!
650        ZQ=ZT
651        QZ0=QS
652!----------------------------------------------------------------------
653        dthv=thvlow-thz0*(0.608*qz0+1.)
654!
655        DU2=MAX((ULOW-UZ0)**2+(VLOW-VZ0)**2,EPSU2)
656        RIB=BTGX*DTHV*ZSL/DU2
657!----------------------------------------------------------------------
658!       IF(RIB>=RIC)THEN
659!         AKMS=MAX( VISC*RDZ,CXCHL)
660!         AKHS=MAX(TVISC*RDZ,CXCHL)
661!----------------------------------------------------------------------
662!       ELSE  !  Turbulent branch
663!----------------------------------------------------------------------
664          ZSLU=ZSL+ZU
665!
666          RZSU=ZSLU/ZU
667!
668          RLOGU=LOG(RZSU)
669
670          ZSLT=ZSL+ZU ! u,v and t are at the same level
671
672          IF ( (IZ0TLND==0) .or. (VEGTYP == ISURBAN) ) THEN
673             ! Just use the original CZIL value.
674             CZIL = 0.1
675          ELSE
676             ! Modify CZIL according to Chen & Zhang, 2009
677             ! CZIL = 10 ** -0.40 H, ( where H = 10*Zo )
678             CZIL = 10.0 ** ( -0.40 * ( Z0 / 0.07 ) )
679          ENDIF
680          ZILFC=-CZIL*VKARMAN*SQVISC
681         
682!----------------------------------------------------------------------
683!
684!mp   Topo modification of ZILFC term
685!       
686!         TOPOTERM=TOPOFAC*ZSFC**2.
687!         TOPOTERM=MAX(TOPOTERM,3.0)
688!
689!  RIB modification to ZILFC term
690!
691!         CZETMAX = 20.
692          CZETMAX = 10.
693! stable
694          IF(DTHV>0.)THEN
695!            ZZIL=ZILFC*TOPOTERM
696            IF (RIB<RIC) THEN
697              ZZIL=ZILFC*(1.0+(RIB/RIC)*(RIB/RIC)*CZETMAX)
698            ELSE
699              ZZIL=ZILFC*(1.0+CZETMAX)
700            ENDIF
701! unstable
702          ELSE
703            ZZIL=ZILFC
704          ENDIF
705!
706!----------------------------------------------------------------------
707!
708          land_point_iteration: DO ITR=1,ITRMX
709!
710!----------------------------------------------------------------------
711!***  ZILITINKEVITCH FIX FOR ZT
712!----------------------------------------------------------------------
713!
714! oldform   ZT=MAX(EXP(ZZIL*SQRT(USTAR*ZU))*ZU,EPSZT)
715            ZT=MAX(EXP(ZZIL*SQRT(USTAR*Z0BASE))*Z0BASE,EPSZT)
716!
717            RZST=ZSLT/ZT
718            RLOGT=LOG(RZST)
719!
720!----------------------------------------------------------------------
721!***  1./MONIN-OBUKHOV LENGTH-SCALE
722!----------------------------------------------------------------------
723!
724            RLMO=ELFC*AKHS*DTHV/USTAR**3
725            ZETALU=ZSLU*RLMO
726            ZETALT=ZSLT*RLMO
727            ZETAU=ZU*RLMO
728            ZETAT=ZT*RLMO
729!
730            ZETALU=MIN(MAX(ZETALU,ZTMIN2),ZTMAX2)
731            ZETALT=MIN(MAX(ZETALT,ZTMIN2),ZTMAX2)
732            ZETAU=MIN(MAX(ZETAU,ZTMIN2/RZSU),ZTMAX2/RZSU)
733            ZETAT=MIN(MAX(ZETAT,ZTMIN2/RZST),ZTMAX2/RZST)
734!
735!----------------------------------------------------------------------
736!***  LAND FUNCTIONS
737!----------------------------------------------------------------------
738!
739            RZ=(ZETAU-ZTMIN2)/DZETA2
740            K=INT(RZ)
741            RDZT=RZ-REAL(K)
742            K=MIN(K,KZTM2)
743            K=MAX(K,0)
744            PSMZ=(PSIM2(K+2)-PSIM2(K+1))*RDZT+PSIM2(K+1)
745!
746            RZ=(ZETALU-ZTMIN2)/DZETA2
747            K=INT(RZ)
748            RDZT=RZ-REAL(K)
749            K=MIN(K,KZTM2)
750            K=MAX(K,0)
751            PSMZL=(PSIM2(K+2)-PSIM2(K+1))*RDZT+PSIM2(K+1)
752!
753            SIMM=PSMZL-PSMZ+RLOGU
754!
755            RZ=(ZETAT-ZTMIN2)/DZETA2
756            K=INT(RZ)
757            RDZT=RZ-REAL(K)
758            K=MIN(K,KZTM2)
759            K=MAX(K,0)
760            PSHZ=(PSIH2(K+2)-PSIH2(K+1))*RDZT+PSIH2(K+1)
761!
762            RZ=(ZETALT-ZTMIN2)/DZETA2
763            K=INT(RZ)
764            RDZT=RZ-REAL(K)
765            K=MIN(K,KZTM2)
766            K=MAX(K,0)
767            PSHZL=(PSIH2(K+2)-PSIH2(K+1))*RDZT+PSIH2(K+1)
768!
769            SIMH=(PSHZL-PSHZ+RLOGT)*FH02
770!----------------------------------------------------------------------
771            USTARK=USTAR*VKARMAN
772            AKMS=MAX(USTARK/SIMM,CXCHL)
773            AKHS=MAX(USTARK/SIMH,CXCHL)
774!
775!----------------------------------------------------------------------
776!***  BELJAARS CORRECTION FOR USTAR
777!----------------------------------------------------------------------
778!
779            IF(DTHV<=0.)THEN                                           !zj
780              WSTAR2=WWST2*ABS(BTGH*AKHS*DTHV)**(2./3.)                !zj
781            ELSE                                                       !zj
782              WSTAR2=0.                                                !zj
783            ENDIF                                                      !zj
784!
785            USTAR=MAX(SQRT(AKMS*SQRT(DU2+WSTAR2)),EPSUST)
786!
787!----------------------------------------------------------------------
788          ENDDO land_point_iteration
789!----------------------------------------------------------------------
790!
791!       ENDIF  !  End of turbulant branch over land
792!
793!----------------------------------------------------------------------
794!
795      ENDIF  !  End of land/sea branch
796!
797!----------------------------------------------------------------------
798!***  COUNTERGRADIENT FIX
799!----------------------------------------------------------------------
800!
801!     HV=-AKHS*DTHV
802!     IF(HV>0.)THEN
803!       FCT=-10.*(BTGX)**(-1./3.)
804!       CT=FCT*(HV/(PBLH*PBLH))**(2./3.)
805!     ELSE
806       CT=0.
807!     ENDIF
808!
809!----------------------------------------------------------------------
810!----------------------------------------------------------------------
811!***  THE FOLLOWING DIAGNOSTIC BLOCK PRODUCES 2-m and 10-m VALUES
812!***  FOR TEMPERATURE, MOISTURE, AND WINDS.  IT IS DONE HERE SINCE
813!***  THE VARIOUS QUANTITIES NEEDED FOR THE COMPUTATION ARE LOST
814!***  UPON EXIT FROM THE ROTUINE.
815!----------------------------------------------------------------------
816!----------------------------------------------------------------------
817!
818      WSTAR=SQRT(WSTAR2)/WWST
819!
820      UMFLX=AKMS*(ULOW -UZ0 )
821      VMFLX=AKMS*(VLOW -VZ0 )
822      HSFLX=AKHS*(THLOW-THZ0)
823      HLFLX=AKHS*(QLOW -QZ0 )
824!----------------------------------------------------------------------
825!     IF(RIB>=RIC)THEN
826!----------------------------------------------------------------------
827!       IF(SEAMASK>0.5)THEN
828!         AKMS10=MAX( VISC/10.,CXCHS)
829!         AKHS02=MAX(TVISC/02.,CXCHS)
830!         AKHS10=MAX(TVISC/10.,CXCHS)
831!       ELSE
832!         AKMS10=MAX( VISC/10.,CXCHL)
833!         AKHS02=MAX(TVISC/02.,CXCHL)
834!         AKHS10=MAX(TVISC/10.,CXCHL)
835!       ENDIF
836!----------------------------------------------------------------------
837!     ELSE
838!----------------------------------------------------------------------
839        ZU10=ZU+10.
840        ZT02=ZT+02.
841        ZT10=ZT+10.
842!
843        RLNU10=LOG(ZU10/ZU)
844        RLNT02=LOG(ZT02/ZT)
845        RLNT10=LOG(ZT10/ZT)
846!
847        ZTAU10=ZU10*RLMO
848        ZTAT02=ZT02*RLMO
849        ZTAT10=ZT10*RLMO
850!
851!----------------------------------------------------------------------
852!***  SEA
853!----------------------------------------------------------------------
854!
855        IF(SEAMASK>0.5)THEN
856!
857!----------------------------------------------------------------------
858          ZTAU10=MIN(MAX(ZTAU10,ZTMIN1),ZTMAX1)
859          ZTAT02=MIN(MAX(ZTAT02,ZTMIN1),ZTMAX1)
860          ZTAT10=MIN(MAX(ZTAT10,ZTMIN1),ZTMAX1)
861!----------------------------------------------------------------------
862          RZ=(ZTAU10-ZTMIN1)/DZETA1
863          K=INT(RZ)
864          RDZT=RZ-REAL(K)
865          K=MIN(K,KZTM2)
866          K=MAX(K,0)
867          PSM10=(PSIM1(K+2)-PSIM1(K+1))*RDZT+PSIM1(K+1)
868!
869          SIMM10=PSM10-PSMZ+RLNU10
870!
871          RZ=(ZTAT02-ZTMIN1)/DZETA1
872          K=INT(RZ)
873          RDZT=RZ-REAL(K)
874          K=MIN(K,KZTM2)
875          K=MAX(K,0)
876          PSH02=(PSIH1(K+2)-PSIH1(K+1))*RDZT+PSIH1(K+1)
877!
878          SIMH02=(PSH02-PSHZ+RLNT02)*FH01
879!
880          RZ=(ZTAT10-ZTMIN1)/DZETA1
881          K=INT(RZ)
882          RDZT=RZ-REAL(K)
883          K=MIN(K,KZTM2)
884          K=MAX(K,0)
885          PSH10=(PSIH1(K+2)-PSIH1(K+1))*RDZT+PSIH1(K+1)
886!
887          SIMH10=(PSH10-PSHZ+RLNT10)*FH01
888!
889          AKMS10=MAX(USTARK/SIMM10,CXCHS)
890          AKHS02=MAX(USTARK/SIMH02,CXCHS)
891          AKHS10=MAX(USTARK/SIMH10,CXCHS)
892!
893!----------------------------------------------------------------------
894!***  LAND
895!----------------------------------------------------------------------
896!
897        ELSE
898!
899!----------------------------------------------------------------------
900          ZTAU10=MIN(MAX(ZTAU10,ZTMIN2),ZTMAX2)
901          ZTAT02=MIN(MAX(ZTAT02,ZTMIN2),ZTMAX2)
902          ZTAT10=MIN(MAX(ZTAT10,ZTMIN2),ZTMAX2)
903!----------------------------------------------------------------------
904          RZ=(ZTAU10-ZTMIN2)/DZETA2
905          K=INT(RZ)
906          RDZT=RZ-REAL(K)
907          K=MIN(K,KZTM2)
908          K=MAX(K,0)
909          PSM10=(PSIM2(K+2)-PSIM2(K+1))*RDZT+PSIM2(K+1)
910!
911          SIMM10=PSM10-PSMZ+RLNU10
912!
913          RZ=(ZTAT02-ZTMIN2)/DZETA2
914          K=INT(RZ)
915          RDZT=RZ-REAL(K)
916          K=MIN(K,KZTM2)
917          K=MAX(K,0)
918          PSH02=(PSIH2(K+2)-PSIH2(K+1))*RDZT+PSIH2(K+1)
919!
920          SIMH02=(PSH02-PSHZ+RLNT02)*FH02
921!
922          RZ=(ZTAT10-ZTMIN2)/DZETA2
923          K=INT(RZ)
924          RDZT=RZ-REAL(K)
925          K=MIN(K,KZTM2)
926          K=MAX(K,0)
927          PSH10=(PSIH2(K+2)-PSIH2(K+1))*RDZT+PSIH2(K+1)
928!
929          SIMH10=(PSH10-PSHZ+RLNT10)*FH02
930!
931          AKMS10=MAX(USTARK/SIMM10,CXCHL)
932          AKHS02=MAX(USTARK/SIMH02,CXCHL)
933          AKHS10=MAX(USTARK/SIMH10,CXCHL)
934!----------------------------------------------------------------------
935        ENDIF
936!----------------------------------------------------------------------
937!     ENDIF
938!----------------------------------------------------------------------
939      U10 =UMFLX/AKMS10+UZ0
940      V10 =VMFLX/AKMS10+VZ0
941      TH02=HSFLX/AKHS02+THZ0
942!
943!***  BE CERTAIN THAT THE 2-M THETA AND 10-M THETA ARE BRACKETED BY
944!***  THE VALUES OF THZ0 AND THLOW.
945!
946      IF(THLOW>THZ0.AND.(TH02<THZ0.OR.TH02>THLOW).OR.                  &
947         THLOW<THZ0.AND.(TH02>THZ0.OR.TH02<THLOW))THEN
948           TH02=THZ0+2.*RDZ*(THLOW-THZ0)
949      ENDIF
950!
951      TH10=HSFLX/AKHS10+THZ0
952!
953      IF(THLOW>THZ0.AND.(TH10<THZ0.OR.TH10>THLOW).OR.                  &
954         THLOW<THZ0.AND.(TH10>THZ0.OR.TH10<THLOW))THEN
955           TH10=THZ0+10.*RDZ*(THLOW-THZ0)
956      ENDIF
957!
958      Q02 =HLFLX/AKHS02+QZ0
959      Q10 =HLFLX/AKHS10+QZ0
960      TERM1=-0.068283/TLOW
961      PSHLTR=PSFC*EXP(TERM1)
962!
963!----------------------------------------------------------------------
964!***  COMPUTE "EQUIVALENT" Z0 TO APPROXIMATE LOCAL SHELTER READINGS.
965!----------------------------------------------------------------------
966!
967      U10E=U10
968      V10E=V10
969!
970      IF(SEAMASK<0.5)THEN
971
972!1st        ZUUZ=MIN(0.5*ZU,0.1)
973!1st        ZU=MAX(0.1*ZU,ZUUZ)
974!tst        ZUUZ=amin1(ZU*0.50,0.3)
975!tst        ZU=amax1(ZU*0.3,ZUUZ)
976
977        ZUUZ=AMIN1(ZU*0.50,0.18)
978        ZU=AMAX1(ZU*0.35,ZUUZ)
979!
980        ZU10=ZU+10.
981        RZSU=ZU10/ZU
982        RLNU10=LOG(RZSU)
983
984        ZETAU=ZU*RLMO
985        ZTAU10=ZU10*RLMO
986
987        ZTAU10=MIN(MAX(ZTAU10,ZTMIN2),ZTMAX2)
988        ZETAU=MIN(MAX(ZETAU,ZTMIN2/RZSU),ZTMAX2/RZSU)
989
990        RZ=(ZTAU10-ZTMIN2)/DZETA2
991        K=INT(RZ)
992        RDZT=RZ-REAL(K)
993        K=MIN(K,KZTM2)
994        K=MAX(K,0)
995        PSM10=(PSIM2(K+2)-PSIM2(K+1))*RDZT+PSIM2(K+1)
996        SIMM10=PSM10-PSMZ+RLNU10
997        EKMS10=MAX(USTARK/SIMM10,CXCHL)
998
999        U10E=UMFLX/EKMS10+UZ0
1000        V10E=VMFLX/EKMS10+VZ0
1001
1002      ENDIF
1003!
1004      U10=U10E
1005      V10=V10E
1006!
1007!----------------------------------------------------------------------
1008!***  SET OTHER WRF DRIVER ARRAYS
1009!----------------------------------------------------------------------
1010!
1011      RLOW=PLOW/(R_D*TLOW)
1012      CHS=AKHS
1013      CHS2=AKHS02
1014      CQS2=AKHS02
1015      HFX=-RLOW*CP*HSFLX
1016      QFX=-RLOW*HLFLX*WETM
1017      FLX_LH=XLV*QFX
1018      FLHC=RLOW*CP*AKHS
1019      FLQC=RLOW*AKHS*WETM
1020!!!   QGH=PQ0/PSHLTR*EXP(A2S*(TSK-A3S)/(TSK-A4S))
1021      QGH=((1.-SEAMASK)*PQ0+SEAMASK*PQ0SEA)                            &
1022     &     /PLOW*EXP(A2S*(TLOW-A3S)/(TLOW-A4S))
1023      QGH=QGH/(1.-QGH)    !Convert to mixing ratio
1024      CPM=CP*(1.+0.8*QLOW)
1025!
1026!***  DO NOT COMPUTE QS OVER LAND POINTS HERE SINCE IT IS
1027!***  A PROGNOSTIC VARIABLE THERE.  IT IS OKAY TO USE IT
1028!***  AS A DIAGNOSTIC OVER WATER SINCE IT WILL CAUSE NO
1029!***  INTERFERENCE BEFORE BEING RECOMPUTED IN MYJPBL.
1030!
1031      IF(SEAMASK>0.5)THEN
1032          QS=PQ0SEA/PSFC                                      &
1033     &         *EXP(A2S*(TSK-A3S)/(TSK-A4S))
1034        QS=QS/(1.-QS)
1035      ENDIF
1036!----------------------------------------------------------------------
1037!
1038      END SUBROUTINE SFCDIF
1039!
1040!----------------------------------------------------------------------
1041      SUBROUTINE MYJSFCINIT(LOWLYR,USTAR,Z0                            &
1042     &                     ,SEAMASK,XICE,IVGTYP,RESTART                &
1043     &                     ,ALLOWED_TO_READ                            &
1044     &                     ,IDS,IDE,JDS,JDE,KDS,KDE                    &
1045     &                     ,IMS,IME,JMS,JME,KMS,KME                    &
1046     &                     ,ITS,ITE,JTS,JTE,KTS,KTE)
1047!----------------------------------------------------------------------
1048      IMPLICIT NONE
1049!----------------------------------------------------------------------
1050      LOGICAL,INTENT(IN) :: RESTART,ALLOWED_TO_READ
1051!
1052      INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                    &
1053     &                     ,IMS,IME,JMS,JME,KMS,KME                    &
1054     &                     ,ITS,ITE,JTS,JTE,KTS,KTE
1055!
1056      INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IVGTYP
1057!
1058      INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: LOWLYR
1059!
1060      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: SEAMASK,XICE
1061!
1062      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: USTAR,Z0
1063!
1064      REAL,DIMENSION(0:30) :: VZ0TBL
1065      REAL,DIMENSION(0:30) :: VZ0TBL_24
1066!
1067      INTEGER :: I,IDUM,IRECV,J,JDUM,K,ITF,JTF,KTF,MAXGBL_IVGTYP       &
1068     &,          MAXLOC_IVGTYP,MPI_COMM_COMP
1069!
1070!     INTEGER :: MPI_INTEGER,MPI_MAX
1071!
1072      REAL :: SM,X,ZETA1,ZETA2,ZRNG1,ZRNG2
1073!
1074      REAL :: PIHF=3.1415926/2.,EPS=1.E-6
1075!----------------------------------------------------------------------
1076      VZ0TBL=                                                          &
1077     &  (/0.,                                                          &
1078     &    2.653,0.826,0.563,1.089,0.854,0.856,0.035,0.238,0.065,0.076  &
1079     &   ,0.011,0.035,0.011,0.000,0.000,0.000,0.000,0.000,0.000,0.000  &
1080     &   ,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000/)
1081
1082      VZ0TBL_24= (/0.,                                                 &
1083     &            1.00,  0.07,  0.07,  0.07,  0.07,  0.15,             &
1084     &            0.08,  0.03,  0.05,  0.86,  0.80,  0.85,             &
1085     &            2.65,  1.09,  0.80,  0.001, 0.04,  0.05,             &
1086     &            0.01,  0.04,  0.06,  0.05,  0.03,  0.001,            &
1087     &            0.000, 0.000, 0.000, 0.000, 0.000, 0.000 /)
1088
1089!----------------------------------------------------------------------
1090!
1091      JTF=MIN0(JTE,JDE-1)
1092      KTF=MIN0(KTE,KDE-1)
1093      ITF=MIN0(ITE,IDE-1)
1094!
1095!
1096!***  FOR NOW, ASSUME SIGMA MODE FOR LOWEST MODEL LAYER
1097!
1098      DO J=JTS,JTF
1099      DO I=ITS,ITF
1100        LOWLYR(I,J)=1
1101!       USTAR(I,J)=EPSUST
1102      ENDDO
1103      ENDDO
1104!----------------------------------------------------------------------
1105#if (NMM_CORE == 1)
1106!
1107      IF(.NOT.RESTART)THEN
1108!       CALL WRF_GET_DM_COMMUNICATOR(MPI_COMM_COMP)
1109!       MAXLOC_IVGTYP=MAXVAL(IVGTYP)
1110!       CALL MPI_ALLREDUCE(MAXLOC_IVGTYP,MAXGBL_IVGTYP,1,MPI_INTEGER   &
1111!    &,                    MPI_MAX,MPI_COMM_COMP,IRECV)
1112        MAXGBL_IVGTYP=MAXVAL(IVGTYP)
1113#ifdef DM_PARALLEL
1114        CALL WRF_DM_MAXVAL(MAXGBL_IVGTYP,IDUM,JDUM)
1115#endif
1116!
1117        IF (MAXGBL_IVGTYP<13) THEN
1118          DO J=JTS,JTE
1119          DO I=ITS,ITE
1120            SM=SEAMASK(I,J)-1.
1121            IF(SM+XICE(I,J)<0.5)THEN
1122              Z0(I,J)=VZ0TBL(IVGTYP(I,J))
1123            ENDIF
1124          ENDDO
1125          ENDDO
1126!
1127        ELSE
1128!
1129          DO J=JTS,JTE
1130          DO I=ITS,ITE
1131            SM=SEAMASK(I,J)-1.
1132            IF(SM+XICE(I,J)<0.5)THEN
1133              Z0(I,J)=VZ0TBL_24(IVGTYP(I,J))
1134            ENDIF
1135          ENDDO
1136          ENDDO
1137!
1138        ENDIF ! Vegtype check
1139!
1140      ENDIF ! Restart check
1141
1142#endif
1143!----------------------------------------------------------------------
1144      IF(.NOT.RESTART)THEN
1145        DO J=JTS,JTE
1146        DO I=ITS,ITF
1147          USTAR(I,J)=0.1
1148        ENDDO
1149        ENDDO
1150      ENDIF
1151!----------------------------------------------------------------------
1152!
1153!***  COMPUTE SURFACE LAYER INTEGRAL FUNCTIONS
1154!
1155!----------------------------------------------------------------------
1156      FH01=1.
1157      FH02=1.
1158!
1159!      ZTMIN1=-10.0
1160!      ZTMAX1=2.0
1161!      ZTMIN2=-10.0
1162!      ZTMAX2=2.0
1163      ZTMIN1=-5.0
1164      ZTMAX1=1.0
1165      ZTMIN2=-5.0
1166      ZTMAX2=1.0
1167!
1168      ZRNG1=ZTMAX1-ZTMIN1
1169      ZRNG2=ZTMAX2-ZTMIN2
1170!
1171      DZETA1=ZRNG1/(KZTM-1)
1172      DZETA2=ZRNG2/(KZTM-1)
1173!
1174!----------------------------------------------------------------------
1175!***  FUNCTION DEFINITION LOOP
1176!----------------------------------------------------------------------
1177!
1178      ZETA1=ZTMIN1
1179      ZETA2=ZTMIN2
1180!
1181      DO K=1,KZTM
1182!
1183!----------------------------------------------------------------------
1184!***  UNSTABLE RANGE
1185!----------------------------------------------------------------------
1186!
1187        IF(ZETA1<0.)THEN
1188!
1189!----------------------------------------------------------------------
1190!***  PAULSON 1970 FUNCTIONS
1191!----------------------------------------------------------------------
1192          X=SQRT(SQRT(1.-16.*ZETA1))
1193!
1194          PSIM1(K)=-2.*LOG((X+1.)/2.)-LOG((X*X+1.)/2.)+2.*ATAN(X)-PIHF
1195          PSIH1(K)=-2.*LOG((X*X+1.)/2.)
1196!
1197!----------------------------------------------------------------------
1198!***  STABLE RANGE
1199!----------------------------------------------------------------------
1200!
1201        ELSE
1202!
1203!----------------------------------------------------------------------
1204!***  PAULSON 1970 FUNCTIONS
1205!----------------------------------------------------------------------
1206!
1207!         PSIM1(K)=5.*ZETA1
1208!         PSIH1(K)=5.*ZETA1
1209!----------------------------------------------------------------------
1210!***   HOLTSLAG AND DE BRUIN 1988
1211!----------------------------------------------------------------------
1212!
1213          PSIM1(K)=0.7*ZETA1+0.75*ZETA1*(6.-0.35*ZETA1)*EXP(-0.35*ZETA1)
1214          PSIH1(K)=0.7*ZETA1+0.75*ZETA1*(6.-0.35*ZETA1)*EXP(-0.35*ZETA1)
1215!----------------------------------------------------------------------
1216!
1217        ENDIF
1218!
1219!----------------------------------------------------------------------
1220!***  UNSTABLE RANGE
1221!----------------------------------------------------------------------
1222!
1223        IF(ZETA2<0.)THEN
1224!
1225!----------------------------------------------------------------------
1226!***  PAULSON 1970 FUNCTIONS
1227!----------------------------------------------------------------------
1228!
1229          X=SQRT(SQRT(1.-16.*ZETA2))
1230!
1231          PSIM2(K)=-2.*LOG((X+1.)/2.)-LOG((X*X+1.)/2.)+2.*ATAN(X)-PIHF
1232          PSIH2(K)=-2.*LOG((X*X+1.)/2.)
1233!----------------------------------------------------------------------
1234!***  STABLE RANGE
1235!----------------------------------------------------------------------
1236!
1237        ELSE
1238!
1239!----------------------------------------------------------------------
1240!***  PAULSON 1970 FUNCTIONS
1241!----------------------------------------------------------------------
1242!
1243!         PSIM2(K)=5.*ZETA2
1244!         PSIH2(K)=5.*ZETA2
1245!
1246!----------------------------------------------------------------------
1247!***  HOLTSLAG AND DE BRUIN 1988
1248!----------------------------------------------------------------------
1249!
1250          PSIM2(K)=0.7*ZETA2+0.75*ZETA2*(6.-0.35*ZETA2)*EXP(-0.35*ZETA2)
1251          PSIH2(K)=0.7*ZETA2+0.75*ZETA2*(6.-0.35*ZETA2)*EXP(-0.35*ZETA2)
1252!----------------------------------------------------------------------
1253!
1254        ENDIF
1255!
1256!----------------------------------------------------------------------
1257        IF(K==KZTM)THEN
1258          ZTMAX1=ZETA1
1259          ZTMAX2=ZETA2
1260        ENDIF
1261!
1262        ZETA1=ZETA1+DZETA1
1263        ZETA2=ZETA2+DZETA2
1264!----------------------------------------------------------------------
1265      ENDDO
1266!----------------------------------------------------------------------
1267      ZTMAX1=ZTMAX1-EPS
1268      ZTMAX2=ZTMAX2-EPS
1269!----------------------------------------------------------------------
1270!
1271      END SUBROUTINE MYJSFCINIT
1272!
1273!----------------------------------------------------------------------
1274!
1275      END MODULE MODULE_SF_MYJSFC
1276!
1277!----------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.