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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

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