source: trunk/WRF.COMMON/WRFV3/phys/module_sf_sfclay.F @ 3431

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

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

File size: 45.1 KB
Line 
1!WRF:MODEL_LAYER:PHYSICS
2!
3MODULE module_sf_sfclay
4
5 REAL    , PARAMETER ::  VCONVC=1.
6 REAL    , PARAMETER ::  CZO=0.0185
7 REAL    , PARAMETER ::  OZO=1.59E-5
8
9 REAL,   DIMENSION(0:1000 ),SAVE          :: PSIMTB,PSIHTB
10
11CONTAINS
12
13!-------------------------------------------------------------------
14   SUBROUTINE SFCLAY(U3D,V3D,T3D,QV3D,P3D,dz8w,                    &
15                     CP,G,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM,      &
16                     ZNT,UST,PBLH,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
17                     XLAND,HFX,QFX,LH,TSK,FLHC,FLQC,QGH,QSFC,RMOL, &
18                     U10,V10,TH2,T2,Q2,                            &
19                     GZ1OZ0,WSPD,BR,ISFFLX,DX,                     &
20                     SVP1,SVP2,SVP3,SVPT0,EP1,EP2,                 &
21                     KARMAN,EOMEG,STBOLT,                          &
22                     P1000mb,                                      &
23                     ids,ide, jds,jde, kds,kde,                    &
24                     ims,ime, jms,jme, kms,kme,                    &
25                     its,ite, jts,jte, kts,kte,                    &
26                     uratx,vratx,tratx,                            &
27                     ustm,ck,cka,cd,cda,isftcflx                   )
28!-------------------------------------------------------------------
29      IMPLICIT NONE
30!-------------------------------------------------------------------
31!-- U3D         3D u-velocity interpolated to theta points (m/s)
32!-- V3D         3D v-velocity interpolated to theta points (m/s)
33!-- T3D         temperature (K)
34!-- QV3D        3D water vapor mixing ratio (Kg/Kg)
35!-- P3D         3D pressure (Pa)
36!-- dz8w        dz between full levels (m)
37!-- CP          heat capacity at constant pressure for dry air (J/kg/K)
38!-- G           acceleration due to gravity (m/s^2)
39!-- ROVCP       R/CP
40!-- R           gas constant for dry air (J/kg/K)
41!-- XLV         latent heat of vaporization for water (J/kg)
42!-- PSFC        surface pressure (Pa)
43!-- ZNT         roughness length (m)
44!-- UST         u* in similarity theory (m/s)
45!-- USTM        u* in similarity theory (m/s) without vconv correction
46!               used to couple with TKE scheme
47!-- PBLH        PBL height from previous time (m)
48!-- MAVAIL      surface moisture availability (between 0 and 1)
49!-- ZOL         z/L height over Monin-Obukhov length
50!-- MOL         T* (similarity theory) (K)
51!-- REGIME      flag indicating PBL regime (stable, unstable, etc.)
52!-- PSIM        similarity stability function for momentum
53!-- PSIH        similarity stability function for heat
54!-- XLAND       land mask (1 for land, 2 for water)
55!-- HFX         upward heat flux at the surface (W/m^2)
56!-- QFX         upward moisture flux at the surface (kg/m^2/s)
57!-- LH          net upward latent heat flux at surface (W/m^2)
58!-- TSK         surface temperature (K)
59!-- FLHC        exchange coefficient for heat (W/m^2/K)
60!-- FLQC        exchange coefficient for moisture (kg/m^2/s)
61!-- CHS         heat/moisture exchange coefficient for LSM (m/s)
62!-- QGH         lowest-level saturated mixing ratio
63!-- QSFC        ground saturated mixing ratio
64!-- uratx       ratio of surface U to U10
65!-- vratx       ratio of surface V to V10
66!-- tratx       ratio of surface T to TH2
67!-- U10         diagnostic 10m u wind
68!-- V10         diagnostic 10m v wind
69!-- TH2         diagnostic 2m theta (K)
70!-- T2          diagnostic 2m temperature (K)
71!-- Q2          diagnostic 2m mixing ratio (kg/kg)
72!-- GZ1OZ0      log(z/z0) where z0 is roughness length
73!-- WSPD        wind speed at lowest model level (m/s)
74!-- BR          bulk Richardson number in surface layer
75!-- ISFFLX      isfflx=1 for surface heat and moisture fluxes
76!-- DX          horizontal grid size (m)
77!-- SVP1        constant for saturation vapor pressure (kPa)
78!-- SVP2        constant for saturation vapor pressure (dimensionless)
79!-- SVP3        constant for saturation vapor pressure (K)
80!-- SVPT0       constant for saturation vapor pressure (K)
81!-- EP1         constant for virtual temperature (R_v/R_d - 1) (dimensionless)
82!-- EP2         constant for specific humidity calculation
83!               (R_d/R_v) (dimensionless)
84!-- KARMAN      Von Karman constant
85!-- EOMEG       angular velocity of earth's rotation (rad/s)
86!-- STBOLT      Stefan-Boltzmann constant (W/m^2/K^4)
87!-- ck          enthalpy exchange coeff at 10 meters
88!-- cd          momentum exchange coeff at 10 meters
89!-- cka         enthalpy exchange coeff at the lowest model level
90!-- cda         momentum exchange coeff at the lowest model level
91!-- isftcflx    =0, (Charnock and Carlson-Boland); =1, AHW Ck, Cd
92!-- ids         start index for i in domain
93!-- ide         end index for i in domain
94!-- jds         start index for j in domain
95!-- jde         end index for j in domain
96!-- kds         start index for k in domain
97!-- kde         end index for k in domain
98!-- ims         start index for i in memory
99!-- ime         end index for i in memory
100!-- jms         start index for j in memory
101!-- jme         end index for j in memory
102!-- kms         start index for k in memory
103!-- kme         end index for k in memory
104!-- its         start index for i in tile
105!-- ite         end index for i in tile
106!-- jts         start index for j in tile
107!-- jte         end index for j in tile
108!-- kts         start index for k in tile
109!-- kte         end index for k in tile
110!-------------------------------------------------------------------
111      INTEGER,  INTENT(IN )   ::        ids,ide, jds,jde, kds,kde, &
112                                        ims,ime, jms,jme, kms,kme, &
113                                        its,ite, jts,jte, kts,kte
114!                                                               
115      INTEGER,  INTENT(IN )   ::        ISFFLX
116      REAL,     INTENT(IN )   ::        SVP1,SVP2,SVP3,SVPT0
117      REAL,     INTENT(IN )   ::        EP1,EP2,KARMAN,EOMEG,STBOLT
118      REAL,     INTENT(IN )   ::        P1000mb
119!
120      REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
121                INTENT(IN   )   ::                           dz8w
122                                       
123      REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
124                INTENT(IN   )   ::                           QV3D, &
125                                                              P3D, &
126                                                              T3D
127
128      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
129                INTENT(IN   )               ::             MAVAIL, &
130                                                             PBLH, &
131                                                            XLAND, &
132                                                              TSK
133      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
134                INTENT(OUT  )               ::                U10, &
135                                                              V10, &
136                                                              TH2, &
137                                                               T2, &
138                                                               Q2, &
139                                                             QSFC
140
141!
142      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
143                INTENT(INOUT)               ::             REGIME, &
144                                                              HFX, &
145                                                              QFX, &
146                                                               LH, &
147                                                          MOL,RMOL
148!m the following 5 are change to memory size
149!
150      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
151                INTENT(INOUT)   ::                 GZ1OZ0,WSPD,BR, &
152                                                        PSIM,PSIH
153
154      REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
155                INTENT(IN   )   ::                            U3D, &
156                                                              V3D
157                                       
158      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
159                INTENT(IN   )               ::               PSFC
160
161      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
162                INTENT(INOUT)   ::                            ZNT, &
163                                                              ZOL, &
164                                                              UST, &
165                                                              CPM, &
166                                                             CHS2, &
167                                                             CQS2, &
168                                                              CHS
169
170      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
171                INTENT(INOUT)   ::                      FLHC,FLQC
172
173      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
174                INTENT(INOUT)   ::                                 &
175                                                              QGH
176
177
178                                   
179      REAL,     INTENT(IN   )               ::   CP,G,ROVCP,R,XLV,DX
180 
181      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme )              , &
182                INTENT(OUT)     ::              uratx,vratx,tratx
183
184      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme )              , &
185                INTENT(OUT)     ::              ck,cka,cd,cda,ustm
186
187      INTEGER,  OPTIONAL,  INTENT(IN )   ::     ISFTCFLX
188
189! LOCAL VARS
190
191      REAL,     DIMENSION( its:ite ) ::                       U1D, &
192                                                              V1D, &
193                                                             QV1D, &
194                                                              P1D, &
195                                                              T1D
196
197      REAL,     DIMENSION( its:ite ) ::                    dz8w1d
198
199      INTEGER ::  I,J
200
201      DO J=jts,jte
202        DO i=its,ite
203          dz8w1d(I) = dz8w(i,1,j)
204        ENDDO
205   
206        DO i=its,ite
207           U1D(i) =U3D(i,1,j)
208           V1D(i) =V3D(i,1,j)
209           QV1D(i)=QV3D(i,1,j)
210           P1D(i) =P3D(i,1,j)
211           T1D(i) =T3D(i,1,j)
212        ENDDO
213
214        !  Sending array starting locations of optional variables may cause
215        !  troubles, so we explicitly change the call.
216
217        CALL SFCLAY1D(J,U1D,V1D,T1D,QV1D,P1D,dz8w1d,               &
218                CP,G,ROVCP,R,XLV,PSFC(ims,j),CHS(ims,j),CHS2(ims,j),&
219                CQS2(ims,j),CPM(ims,j),PBLH(ims,j), RMOL(ims,j),   &
220                ZNT(ims,j),UST(ims,j),MAVAIL(ims,j),ZOL(ims,j),    &
221                MOL(ims,j),REGIME(ims,j),PSIM(ims,j),PSIH(ims,j),  &
222                XLAND(ims,j),HFX(ims,j),QFX(ims,j),TSK(ims,j),     &
223                U10(ims,j),V10(ims,j),TH2(ims,j),T2(ims,j),        &
224                Q2(ims,j),FLHC(ims,j),FLQC(ims,j),QGH(ims,j),      &
225                QSFC(ims,j),LH(ims,j),                             &
226                GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),ISFFLX,DX,     &
227                SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,  &
228                P1000mb,                                           &
229                ids,ide, jds,jde, kds,kde,                         &
230                ims,ime, jms,jme, kms,kme,                         &
231                its,ite, jts,jte, kts,kte                          &
232#if ( EM_CORE == 1 )
233                ,uratx(ims,j),vratx(ims,j),tratx(ims,j),           &
234                isftcflx,                                          &
235                USTM(ims,j),CK(ims,j),CKA(ims,j),                  &
236                CD(ims,j),CDA(ims,j)                               &
237#endif
238                                                                   )
239      ENDDO
240
241
242   END SUBROUTINE SFCLAY
243
244
245!-------------------------------------------------------------------
246   SUBROUTINE SFCLAY1D(J,UX,VX,T1D,QV1D,P1D,dz8w1d,                &
247                     CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2,CPM,PBLH,RMOL, &
248                     ZNT,UST,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH,      &
249                     XLAND,HFX,QFX,TSK,                            &
250                     U10,V10,TH2,T2,Q2,FLHC,FLQC,QGH,              &
251                     QSFC,LH,GZ1OZ0,WSPD,BR,ISFFLX,DX,             &
252                     SVP1,SVP2,SVP3,SVPT0,EP1,EP2,                 &
253                     KARMAN,EOMEG,STBOLT,                          &
254                     P1000mb,                                      &
255                     ids,ide, jds,jde, kds,kde,                    &
256                     ims,ime, jms,jme, kms,kme,                    &
257                     its,ite, jts,jte, kts,kte,                    &
258                     uratx,vratx,tratx,                            &
259                     isftcflx,                                     &
260                     ustm,ck,cka,cd,cda                            )
261!-------------------------------------------------------------------
262      IMPLICIT NONE
263!-------------------------------------------------------------------
264      REAL,     PARAMETER     ::        XKA=2.4E-5
265      REAL,     PARAMETER     ::        PRT=1.
266
267      INTEGER,  INTENT(IN )   ::        ids,ide, jds,jde, kds,kde, &
268                                        ims,ime, jms,jme, kms,kme, &
269                                        its,ite, jts,jte, kts,kte, &
270                                        J
271!                                                               
272      INTEGER,  INTENT(IN )   ::        ISFFLX
273      REAL,     INTENT(IN )   ::        SVP1,SVP2,SVP3,SVPT0
274      REAL,     INTENT(IN )   ::        EP1,EP2,KARMAN,EOMEG,STBOLT
275      REAL,     INTENT(IN )   ::        P1000mb
276
277!
278      REAL,     DIMENSION( ims:ime )                             , &
279                INTENT(IN   )               ::             MAVAIL, &
280                                                             PBLH, &
281                                                            XLAND, &
282                                                              TSK
283!
284      REAL,     DIMENSION( ims:ime )                             , &
285                INTENT(IN   )               ::             PSFCPA
286
287      REAL,     DIMENSION( ims:ime )                             , &
288                INTENT(INOUT)               ::             REGIME, &
289                                                              HFX, &
290                                                              QFX, &
291                                                         MOL,RMOL
292!m the following 5 are changed to memory size---
293!
294      REAL,     DIMENSION( ims:ime )                             , &
295                INTENT(INOUT)   ::                 GZ1OZ0,WSPD,BR, &
296                                                        PSIM,PSIH
297
298      REAL,     DIMENSION( ims:ime )                             , &
299                INTENT(INOUT)   ::                            ZNT, &
300                                                              ZOL, &
301                                                              UST, &
302                                                              CPM, &
303                                                             CHS2, &
304                                                             CQS2, &
305                                                              CHS
306
307      REAL,     DIMENSION( ims:ime )                             , &
308                INTENT(INOUT)   ::                      FLHC,FLQC
309
310      REAL,     DIMENSION( ims:ime )                             , &
311                INTENT(INOUT)   ::                                 &
312                                                              QGH
313
314      REAL,     DIMENSION( ims:ime )                             , &
315                INTENT(OUT)     ::                        U10,V10, &
316                                                TH2,T2,Q2,QSFC,LH
317
318                                   
319      REAL,     INTENT(IN   )               ::   CP,G,ROVCP,R,XLV,DX
320
321! MODULE-LOCAL VARIABLES, DEFINED IN SUBROUTINE SFCLAY
322      REAL,     DIMENSION( its:ite ),  INTENT(IN   )   ::  dz8w1d
323
324      REAL,     DIMENSION( its:ite ),  INTENT(IN   )   ::      UX, &
325                                                               VX, &
326                                                             QV1D, &
327                                                              P1D, &
328                                                              T1D
329 
330      REAL, OPTIONAL, DIMENSION( ims:ime )                       , &
331                INTENT(OUT)     ::              uratx,vratx,tratx
332
333      REAL, OPTIONAL, DIMENSION( ims:ime )                       , &
334                INTENT(OUT)     ::              ck,cka,cd,cda,ustm
335
336      INTEGER,  OPTIONAL,  INTENT(IN )   ::     ISFTCFLX
337
338! LOCAL VARS
339
340      REAL,     DIMENSION( its:ite )        ::                 ZA, &
341                                                        THVX,ZQKL, &
342                                                           ZQKLP1, &
343                                                           THX,QX, &
344                                                            PSIH2, &
345                                                            PSIM2, &
346                                                           PSIH10, &
347                                                           PSIM10, &
348                                                           DENOMQ, &
349                                                            WSPDI, &
350                                                           GZ2OZ0, &
351                                                           GZ10OZ0
352!
353      REAL,     DIMENSION( its:ite )        ::                     &
354                                                      RHOX,GOVRTH, &
355                                                            TGDSA
356!
357      REAL,     DIMENSION( its:ite)         ::          SCR3,SCR4
358      REAL,     DIMENSION( its:ite )        ::         THGB, PSFC
359!
360      INTEGER                               ::                 KL
361
362      INTEGER ::  N,I,K,KK,L,NZOL,NK,NZOL2,NZOL10
363
364      REAL    ::  PL,THCON,TVCON,E1
365      REAL    ::  ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL,RZOL2,RZOL10,ZOL2,ZOL10
366      REAL    ::  DTG,PSIX,DTTHX,PSIX10,PSIT,PSIT2,PSIQ,PSIQ2,PSIQ10
367      REAL    ::  FLUXC,VSGD,Z0Q
368!-------------------------------------------------------------------
369      KL=kte
370
371      DO i=its,ite
372! PSFC cb
373         PSFC(I)=PSFCPA(I)/1000.
374      ENDDO
375!                                                     
376!----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE: 
377!                                                           
378      DO 5 I=its,ite                                   
379        TGDSA(I)=TSK(I)                                   
380! PSFC cb
381!        THGB(I)=TSK(I)*(100./PSFC(I))**ROVCP               
382        THGB(I)=TSK(I)*(P1000mb/PSFCPA(I))**ROVCP   
383    5 CONTINUE                                               
384!                                                           
385!-----DECOUPLE FLUX-FORM VARIABLES TO GIVE U,V,T,THETA,THETA-VIR.,
386!     T-VIR., QV, AND QC AT CROSS POINTS AND AT KTAU-1. 
387!                                                                 
388!     *** NOTE ***                                           
389!         THE BOUNDARY WINDS MAY NOT BE ADEQUATELY AFFECTED BY FRICTION,         
390!         SO USE ONLY INTERIOR VALUES OF UX AND VX TO CALCULATE
391!         TENDENCIES.                             
392!                                                           
393   10 CONTINUE                                                     
394
395!     DO 24 I=its,ite
396!        UX(I)=U1D(I)
397!        VX(I)=V1D(I)
398!  24 CONTINUE                                             
399                                                             
400   26 CONTINUE                                               
401                                                   
402!.....SCR3(I,K) STORE TEMPERATURE,                           
403!     SCR4(I,K) STORE VIRTUAL TEMPERATURE.                                       
404                                                                                 
405      DO 30 I=its,ite
406! PL cb
407         PL=P1D(I)/1000.
408         SCR3(I)=T1D(I)                                                   
409!         THCON=(100./PL)**ROVCP                                                 
410         THCON=(P1000mb*0.001/PL)**ROVCP
411         THX(I)=SCR3(I)*THCON                                               
412         SCR4(I)=SCR3(I)                                                   
413         THVX(I)=THX(I)                                                     
414         QX(I)=0.                                                             
415   30 CONTINUE                                                                 
416!                                                                               
417      DO I=its,ite
418         QGH(I)=0.                                                               
419         FLHC(I)=0.                                                               
420         FLQC(I)=0.                                                               
421         CPM(I)=CP                                                               
422      ENDDO
423!                                                                               
424!     IF(IDRY.EQ.1)GOTO 80                                                   
425      DO 50 I=its,ite
426         QX(I)=QV1D(I)                                                   
427         TVCON=(1.+EP1*QX(I))                                     
428         THVX(I)=THX(I)*TVCON                                               
429         SCR4(I)=SCR3(I)*TVCON                                             
430   50 CONTINUE                                                                 
431!                                                                               
432      DO 60 I=its,ite
433        E1=SVP1*EXP(SVP2*(TGDSA(I)-SVPT0)/(TGDSA(I)-SVP3))                       
434        QSFC(I)=EP2*E1/(PSFC(I)-E1)                                                 
435! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP CONSISTENT WITH MYJSFC CHANGE
436! Q2SAT = QGH IN LSM
437        E1=SVP1*EXP(SVP2*(T1D(I)-SVPT0)/(T1D(I)-SVP3))                       
438        PL=P1D(I)/1000.
439        QGH(I)=EP2*E1/(PL-E1)                                                 
440        CPM(I)=CP*(1.+0.8*QX(I))                                   
441   60 CONTINUE                                                                   
442   80 CONTINUE
443                                                                                 
444!-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND             
445!     LEVEL, AND THE LAYER THICKNESSES.                                         
446                                                                                 
447      DO 90 I=its,ite
448        ZQKLP1(I)=0.
449        RHOX(I)=PSFC(I)*1000./(R*SCR4(I))                                       
450   90 CONTINUE                                                                   
451!                                                                               
452      DO 110 I=its,ite                                                   
453           ZQKL(I)=dz8w1d(I)+ZQKLP1(I)
454  110 CONTINUE                                                                 
455!                                                                               
456      DO 120 I=its,ite
457         ZA(I)=0.5*(ZQKL(I)+ZQKLP1(I))                                       
458  120 CONTINUE                                                                 
459!                                                                               
460      DO 160 I=its,ite
461        GOVRTH(I)=G/THX(I)                                                   
462  160 CONTINUE                                                                   
463                                                                                 
464!-----CALCULATE BULK RICHARDSON NO. OF SURFACE LAYER, ACCORDING TO               
465!     AKB(1976), EQ(12).                                                         
466                   
467      DO 260 I=its,ite
468        GZ1OZ0(I)=ALOG(ZA(I)/ZNT(I))                                       
469        GZ2OZ0(I)=ALOG(2./ZNT(I))                                       
470        GZ10OZ0(I)=ALOG(10./ZNT(I))                                       
471        IF((XLAND(I)-1.5).GE.0)THEN                                           
472          ZL=ZNT(I)                                                           
473        ELSE                                                                     
474          ZL=0.01                                                               
475        ENDIF                                                                   
476        WSPD(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I))                       
477
478        TSKV=THGB(I)*(1.+EP1*QSFC(I)*MAVAIL(I))                     
479        DTHVDZ=(THVX(I)-TSKV)                                                 
480!  Convective velocity scale Vc and subgrid-scale velocity Vsg
481!  following Beljaars (1995, QJRMS) and Mahrt and Sun (1995, MWR)
482!                                ... HONG Aug. 2001
483!
484!       VCONV = 0.25*sqrt(g/tskv*pblh(i)*dthvm)
485!      Use Beljaars over land, old MM5 (Wyngaard) formula over water
486        if (xland(i).lt.1.5) then
487        fluxc = max(hfx(i)/rhox(i)/cp                    &
488              + ep1*tskv*qfx(i)/rhox(i),0.)
489        VCONV = vconvc*(g/tgdsa(i)*pblh(i)*fluxc)**.33
490        else
491        IF(-DTHVDZ.GE.0)THEN
492          DTHVM=-DTHVDZ
493        ELSE
494          DTHVM=0.
495        ENDIF
496        VCONV = 2.*SQRT(DTHVM)
497        endif
498! Mahrt and Sun low-res correction
499        VSGD = 0.32 * (max(dx/5000.-1.,0.))**.33
500        WSPD(I)=SQRT(WSPD(I)*WSPD(I)+VCONV*VCONV+vsgd*vsgd)
501        WSPD(I)=AMAX1(WSPD(I),0.1)
502        BR(I)=GOVRTH(I)*ZA(I)*DTHVDZ/(WSPD(I)*WSPD(I))                       
503!  IF PREVIOUSLY UNSTABLE, DO NOT LET INTO REGIMES 1 AND 2
504        IF(MOL(I).LT.0.)BR(I)=AMIN1(BR(I),0.0)
505!jdf
506        RMOL(I)=-GOVRTH(I)*DTHVDZ*ZA(I)*KARMAN
507!jdf
508
509  260 CONTINUE                                                                   
510
511!                                                                               
512!-----DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS:           
513!                                                                               
514!                                                                               
515!     THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.)           
516!     AND HOL (HEIGHT OF PBL/MONIN-OBUKHOV LENGTH).                             
517!                                                                               
518!     CRITERIA FOR THE CLASSES ARE AS FOLLOWS:                                   
519!                                                                               
520!        1. BR .GE. 0.2;                                                         
521!               REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1),               
522!                                                                               
523!        2. BR .LT. 0.2 .AND. BR .GT. 0.0;                                       
524!               REPRESENTS DAMPED MECHANICAL TURBULENT CONDITIONS               
525!               (REGIME=2),                                                     
526!                                                                               
527!        3. BR .EQ. 0.0                                                         
528!               REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3),             
529!                                                                               
530!        4. BR .LT. 0.0                                                         
531!               REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4).               
532!                                                                               
533!CCCCC                                                                           
534
535      DO 320 I=its,ite
536!CCCCC                                                                           
537!CC     REMOVE REGIME 3 DEPENDENCE ON PBL HEIGHT                                 
538!CC          IF(BR(I).LT.0..AND.HOL(I,J).GT.1.5)GOTO 310                         
539        IF(BR(I).LT.0.)GOTO 310                                                 
540!                                                                               
541!-----CLASS 1; STABLE (NIGHTTIME) CONDITIONS:                                   
542!                                                                               
543        IF(BR(I).LT.0.2)GOTO 270                                                 
544        REGIME(I)=1.                                                           
545        PSIM(I)=-10.*GZ1OZ0(I)                                                   
546!    LOWER LIMIT ON PSI IN STABLE CONDITIONS                                     
547        PSIM(I)=AMAX1(PSIM(I),-10.)                                             
548        PSIH(I)=PSIM(I)                                                         
549        PSIM10(I)=10./ZA(I)*PSIM(I)
550        PSIM10(I)=AMAX1(PSIM10(I),-10.)                               
551        PSIH10(I)=PSIM10(I)                                         
552        PSIM2(I)=2./ZA(I)*PSIM(I)
553        PSIM2(I)=AMAX1(PSIM2(I),-10.)                             
554        PSIH2(I)=PSIM2(I)                                         
555
556!       1.0 over Monin-Obukhov length
557        IF(UST(I).LT.0.01)THEN
558           RMOL(I)=BR(I)*GZ1OZ0(I) !ZA/L
559        ELSE
560           RMOL(I)=KARMAN*GOVRTH(I)*ZA(I)*MOL(I)/(UST(I)*UST(I)) !ZA/L
561        ENDIF
562        RMOL(I)=AMIN1(RMOL(I),9.999) ! ZA/L
563        RMOL(I) = RMOL(I)/ZA(I) !1.0/L
564
565        GOTO 320                                                                 
566!                                                                               
567!-----CLASS 2; DAMPED MECHANICAL TURBULENCE:                                     
568!                                                                               
569  270   IF(BR(I).EQ.0.0)GOTO 280                                                 
570        REGIME(I)=2.                                                           
571        PSIM(I)=-5.0*BR(I)*GZ1OZ0(I)/(1.1-5.0*BR(I))                             
572!    LOWER LIMIT ON PSI IN STABLE CONDITIONS                                     
573        PSIM(I)=AMAX1(PSIM(I),-10.)                                             
574!.....AKB(1976), EQ(16).                                                         
575        PSIH(I)=PSIM(I)                                                         
576        PSIM10(I)=10./ZA(I)*PSIM(I)
577        PSIM10(I)=AMAX1(PSIM10(I),-10.)                               
578        PSIH10(I)=PSIM10(I)                                         
579        PSIM2(I)=2./ZA(I)*PSIM(I)
580        PSIM2(I)=AMAX1(PSIM2(I),-10.)                             
581        PSIH2(I)=PSIM2(I)                                         
582
583        ! Linear form: PSIM = -0.5*ZA/L; e.g, see eqn 16 of
584        ! Blackadar, Modeling the nocturnal boundary layer, Preprints,
585        ! Third Symposium on Atmospheric Turbulence Diffusion and Air Quality,
586        ! Raleigh, NC, 1976
587        ZOL(I) = BR(I)*GZ1OZ0(I)/(1.00001-5.0*BR(I))
588
589        if ( ZOL(I) .GT. 0.5 ) then ! linear form ok
590           ! Holtslag and de Bruin, J. App. Meteor 27, 689-704, 1988;
591           ! see also, Launiainen, Boundary-Layer Meteor 76,165-179, 1995
592           ! Eqn (8) of Launiainen, 1995
593           ZOL(I) = ( 1.89*GZ1OZ0(I) + 44.2 ) * BR(I)*BR(I)    &
594                + ( 1.18*GZ1OZ0(I) - 1.37 ) * BR(I)
595           ZOL(I)=AMIN1(ZOL(I),9.999)
596        end if
597
598        ! 1.0 over Monin-Obukhov length
599        RMOL(I)= ZOL(I)/ZA(I)
600
601        GOTO 320                                                                 
602!                                                                               
603!-----CLASS 3; FORCED CONVECTION:                                               
604!                                                                               
605  280   REGIME(I)=3.                                                           
606        PSIM(I)=0.0                                                             
607        PSIH(I)=PSIM(I)                                                         
608        PSIM10(I)=0.                                                   
609        PSIH10(I)=PSIM10(I)                                           
610        PSIM2(I)=0.                                                 
611        PSIH2(I)=PSIM2(I)                                           
612
613                                                                                 
614        IF(UST(I).LT.0.01)THEN                                                 
615          ZOL(I)=BR(I)*GZ1OZ0(I)                                               
616        ELSE                                                                     
617          ZOL(I)=KARMAN*GOVRTH(I)*ZA(I)*MOL(I)/(UST(I)*UST(I))
618        ENDIF                                                                   
619
620        RMOL(I) = ZOL(I)/ZA(I) 
621
622        GOTO 320                                                                 
623!                                                                               
624!-----CLASS 4; FREE CONVECTION:                                                 
625!                                                                               
626  310   CONTINUE                                                                 
627        REGIME(I)=4.                                                           
628        IF(UST(I).LT.0.01)THEN                                                 
629          ZOL(I)=BR(I)*GZ1OZ0(I)                                               
630        ELSE                                                                     
631          ZOL(I)=KARMAN*GOVRTH(I)*ZA(I)*MOL(I)/(UST(I)*UST(I))
632        ENDIF                                                                   
633        ZOL10=10./ZA(I)*ZOL(I)                                   
634        ZOL2=2./ZA(I)*ZOL(I)                                     
635        ZOL(I)=AMIN1(ZOL(I),0.)                                             
636        ZOL(I)=AMAX1(ZOL(I),-9.9999)                                         
637        ZOL10=AMIN1(ZOL10,0.)                                         
638        ZOL10=AMAX1(ZOL10,-9.9999)                                   
639        ZOL2=AMIN1(ZOL2,0.)                                         
640        ZOL2=AMAX1(ZOL2,-9.9999)                                   
641        NZOL=INT(-ZOL(I)*100.)                                                 
642        RZOL=-ZOL(I)*100.-NZOL                                                 
643        NZOL10=INT(-ZOL10*100.)                                       
644        RZOL10=-ZOL10*100.-NZOL10                                     
645        NZOL2=INT(-ZOL2*100.)                                       
646        RZOL2=-ZOL2*100.-NZOL2                                     
647        PSIM(I)=PSIMTB(NZOL)+RZOL*(PSIMTB(NZOL+1)-PSIMTB(NZOL))                 
648        PSIH(I)=PSIHTB(NZOL)+RZOL*(PSIHTB(NZOL+1)-PSIHTB(NZOL))                 
649        PSIM10(I)=PSIMTB(NZOL10)+RZOL10*(PSIMTB(NZOL10+1)-PSIMTB(NZOL10))                                                   
650        PSIH10(I)=PSIHTB(NZOL10)+RZOL10*(PSIHTB(NZOL10+1)-PSIHTB(NZOL10))
651        PSIM2(I)=PSIMTB(NZOL2)+RZOL2*(PSIMTB(NZOL2+1)-PSIMTB(NZOL2))   
652        PSIH2(I)=PSIHTB(NZOL2)+RZOL2*(PSIHTB(NZOL2+1)-PSIHTB(NZOL2))   
653
654!---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND HIGH ROUGHNESS           
655!---  THIS PREVENTS DENOMINATOR IN FLUXES FROM GETTING TOO SMALL                 
656!       PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))                                     
657!       PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))                                     
658        PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
659        PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
660        PSIH2(I)=AMIN1(PSIH2(I),0.9*GZ2OZ0(I))
661        PSIM10(I)=AMIN1(PSIM10(I),0.9*GZ10OZ0(I))
662! AHW: mods to compute ck, cd
663        PSIH10(I)=AMIN1(PSIH10(I),0.9*GZ10OZ0(I))
664
665        RMOL(I) = ZOL(I)/ZA(I) 
666
667  320 CONTINUE                                                                   
668!                                                                               
669!-----COMPUTE THE FRICTIONAL VELOCITY:                                           
670!     ZA(1982) EQS(2.60),(2.61).                                                 
671!                                                                               
672      DO 330 I=its,ite
673        DTG=THX(I)-THGB(I)                                                   
674        PSIX=GZ1OZ0(I)-PSIM(I)                                                   
675        PSIX10=GZ10OZ0(I)-PSIM10(I)
676!     LOWER LIMIT ADDED TO PREVENT LARGE FLHC IN SOIL MODEL
677!     ACTIVATES IN UNSTABLE CONDITIONS WITH THIN LAYERS OR HIGH Z0
678        PSIT=AMAX1(GZ1OZ0(I)-PSIH(I),2.)
679
680        IF((XLAND(I)-1.5).GE.0)THEN                                           
681          ZL=ZNT(I)                                                           
682        ELSE                                                                     
683          ZL=0.01                                                               
684        ENDIF                                                                   
685        PSIQ=ALOG(KARMAN*UST(I)*ZA(I)/XKA+ZA(I)/ZL)-PSIH(I)   
686        PSIT2=GZ2OZ0(I)-PSIH2(I)                                     
687        PSIQ2=ALOG(KARMAN*UST(I)*2./XKA+2./ZL)-PSIH2(I)                                   
688! AHW: mods to compute ck, cd
689        PSIQ10=ALOG(KARMAN*UST(I)*10./XKA+10./ZL)-PSIH10(I)
690        IF ( PRESENT(ISFTCFLX) ) THEN
691           IF ( ISFTCFLX.EQ.1 .AND. (XLAND(I)-1.5).GE.0. ) THEN
692              Z0Q = 1.e-4 + 1.e-3*(MAX(0.,UST(I)-1.))**2
693!             Z0Q = 0.62*1.5E-5/UST(I) + 1.E-3*(MAX(0.,UST(I)-1.))**2
694              PSIQ=ALOG(ZA(I)/Z0Q)-PSIH(I)
695              PSIT=PSIQ
696              PSIQ2=ALOG(2./Z0Q)-PSIH2(I)
697              PSIQ10=ALOG(10./Z0Q)-PSIH10(I)
698              PSIT2=PSIQ2
699           ENDIF
700        ENDIF
701        IF(PRESENT(ck) .and. PRESENT(cd) .and. PRESENT(cka) .and. PRESENT(cda)) THEN
702           Ck(I)=(karman/psix10)*(karman/psiq10)
703           Cd(I)=(karman/psix10)*(karman/psix10)
704           Cka(I)=(karman/psix)*(karman/psiq)
705           Cda(I)=(karman/psix)*(karman/psix)
706        ENDIF
707! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE
708        UST(I)=0.5*UST(I)+0.5*KARMAN*WSPD(I)/PSIX                                             
709! TKE coupling: compute ust without vconv for use in tke scheme
710        WSPDI(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I))
711        IF ( PRESENT(USTM) ) THEN
712        USTM(I)=0.5*USTM(I)+0.5*KARMAN*WSPDI(I)/PSIX
713        ENDIF
714        U10(I)=UX(I)*PSIX10/PSIX                                   
715        V10(I)=VX(I)*PSIX10/PSIX                                   
716        TH2(I)=THGB(I)+DTG*PSIT2/PSIT                               
717        Q2(I)=QSFC(I)+(QX(I)-QSFC(I))*PSIQ2/PSIQ                   
718!        T2(I) = TH2(I)*(PSFC(I)/100.)**ROVCP                     
719        T2(I) = TH2(I)*(PSFCPA(I)/P1000mb)**ROVCP                     
720!       LATER Q2 WILL BE OVERWRITTEN FOR LAND POINTS IN SURFCE     
721!       QA2(I,J) = Q2(I)                                         
722!       UA10(I,J) = U10(I)                                     
723!       VA10(I,J) = V10(I)                                     
724!       write(*,1002)UST(I),KARMAN*WSPD(I),PSIX,KARMAN*WSPD(I)/PSIX
725!                                                                               
726        IF(PRESENT(uratx) .and. PRESENT(vratx) .and. PRESENT(tratx))THEN
727          IF(ABS(U10(I)) .GT. 1.E-10) THEN
728             uratx(I) = UX(I)/U10(I)
729          ELSE
730             uratx(I) = 1.2
731          END IF
732          IF(ABS(V10(I)) .GT. 1.E-10) THEN
733             vratx(I) = VX(I)/V10(I)
734          ELSE
735             vratx(I) = 1.2
736          END IF
737          tratx(I) = THX(I)/TH2(I)
738        ENDIF
739
740        IF((XLAND(I)-1.5).LT.0.)THEN                                           
741          UST(I)=AMAX1(UST(I),0.1)
742        ENDIF                                                                   
743        MOL(I)=KARMAN*DTG/PSIT/PRT                             
744        DENOMQ(I)=PSIQ
745  330 CONTINUE                                                                   
746!                                                                               
747  335 CONTINUE                                                                   
748                                                                                 
749!-----COMPUTE THE SURFACE SENSIBLE AND LATENT HEAT FLUXES:                       
750                                                                                 
751      DO i=its,ite
752        QFX(i)=0.                                                             
753        HFX(i)=0.                                                             
754      ENDDO
755
756      IF (ISFFLX.EQ.0) GOTO 410                                               
757                                                                                 
758!-----OVER WATER, ALTER ROUGHNESS LENGTH (ZNT) ACCORDING TO WIND (UST).         
759                                                                                 
760      DO 360 I=its,ite
761        IF((XLAND(I)-1.5).GE.0)THEN                                           
762          ZNT(I)=CZO*UST(I)*UST(I)/G+OZO                                   
763! AHW: change roughness length, and hence the drag coefficients Ck and Cd
764          IF ( PRESENT(ISFTCFLX) ) THEN
765             IF ( ISFTCFLX.EQ.1 ) THEN
766                ZNT(I)=10.*exp(-9.*UST(I)**(-.3333))
767                ZNT(I)=MIN(ZNT(I),2.85e-3)
768                ZNT(I)=ZNT(I) + 0.11*1.5E-5/AMAX1(UST(I),0.01)
769!               ZNT(I)=MAX(ZNT(I),1.27e-7)
770             ENDIF
771          ENDIF
772          ZL = ZNT(I)
773        ELSE
774          ZL = 0.01
775        ENDIF                                                                   
776        FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/DENOMQ(I)
777!       FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/(   &
778!               ALOG(KARMAN*UST(I)*ZA(I)/XKA+ZA(I)/ZL)-PSIH(I))
779        DTTHX=ABS(THX(I)-THGB(I))                                           
780        IF(DTTHX.GT.1.E-5)THEN                                                   
781          FLHC(I)=CPM(I)*RHOX(I)*UST(I)*MOL(I)/(THX(I)-THGB(I))         
782!         write(*,1001)FLHC(I),CPM(I),RHOX(I),UST(I),MOL(I),THX(I),THGB(I),I
783 1001   format(f8.5,2x,f12.7,2x,f12.10,2x,f12.10,2x,f13.10,2x,f12.8,f12.8,2x,i3)
784        ELSE                                                                     
785          FLHC(I)=0.                                                             
786        ENDIF                                                                   
787  360 CONTINUE                                                                   
788
789!                                                                               
790!-----COMPUTE SURFACE MOIST FLUX:                                               
791!                                                                               
792!     IF(IDRY.EQ.1)GOTO 390                                               
793!                                                                               
794      DO 370 I=its,ite
795        QFX(I)=FLQC(I)*(QSFC(I)-QX(I))                                     
796        QFX(I)=AMAX1(QFX(I),0.)                                           
797        LH(I)=XLV*QFX(I)
798  370 CONTINUE                                                                 
799                                                                               
800!-----COMPUTE SURFACE HEAT FLUX:                                                 
801!                                                                               
802  390 CONTINUE                                                                 
803      DO 400 I=its,ite
804        IF(XLAND(I)-1.5.GT.0.)THEN                                           
805          HFX(I)=FLHC(I)*(THGB(I)-THX(I))                               
806        ELSEIF(XLAND(I)-1.5.LT.0.)THEN                                       
807          HFX(I)=FLHC(I)*(THGB(I)-THX(I))                               
808          HFX(I)=AMAX1(HFX(I),-250.)                                       
809        ENDIF                                                                 
810  400 CONTINUE                                                                 
811         
812      DO I=its,ite
813         IF((XLAND(I)-1.5).GE.0)THEN
814           ZL=ZNT(I)
815         ELSE
816           ZL=0.01
817         ENDIF
818         CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) &
819                /XKA+ZA(I)/ZL)-PSIH(I))
820!        GZ2OZ0(I)=ALOG(2./ZNT(I))
821!        PSIM2(I)=-10.*GZ2OZ0(I)
822!        PSIM2(I)=AMAX1(PSIM2(I),-10.)
823!        PSIH2(I)=PSIM2(I)
824         CQS2(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*2.0  &
825               /XKA+2.0/ZL)-PSIH2(I))
826         CHS2(I)=UST(I)*KARMAN/(GZ2OZ0(I)-PSIH2(I))
827      ENDDO
828                                                                       
829  410 CONTINUE                                                                   
830!jdf
831!     DO I=its,ite
832!       IF(UST(I).GE.0.1) THEN
833!         RMOL(I)=RMOL(I)*(-FLHC(I))/(UST(I)*UST(I)*UST(I))
834!       ELSE
835!         RMOL(I)=RMOL(I)*(-FLHC(I))/(0.1*0.1*0.1)
836!       ENDIF
837!     ENDDO
838!jdf
839
840!                                                                               
841   END SUBROUTINE SFCLAY1D
842
843!====================================================================
844   SUBROUTINE sfclayinit( allowed_to_read )         
845
846   LOGICAL , INTENT(IN)      ::      allowed_to_read
847   INTEGER                   ::      N
848   REAL                      ::      ZOLN,X,Y
849
850   DO N=0,1000
851      ZOLN=-FLOAT(N)*0.01
852      X=(1-16.*ZOLN)**0.25
853      PSIMTB(N)=2*ALOG(0.5*(1+X))+ALOG(0.5*(1+X*X))- &
854                2.*ATAN(X)+2.*ATAN(1.)
855      Y=(1-16*ZOLN)**0.5
856      PSIHTB(N)=2*ALOG(0.5*(1+Y))
857   ENDDO
858
859   END SUBROUTINE sfclayinit
860
861!-------------------------------------------------------------------         
862
863END MODULE module_sf_sfclay
Note: See TracBrowser for help on using the repository browser.