source: lmdz_wrf/trunk/WRFV3/phys/module_sf_mynn.F @ 354

Last change on this file since 354 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:

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