source: trunk/WRF.COMMON/WRFV3/phys/module_sf_myjsfc.F @ 3094

Last change on this file since 3094 was 2759, checked in by aslmd, 2 years ago

adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size

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