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

Last change on this file since 3094 was 11, checked in by aslmd, 14 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien

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