Newer
Older
END IF
IF( LBUDGET_RC ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_RC), 'HTURB', PRRS(:, :, :, 2) )
IF( LBUDGET_RI ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_RI), 'HTURB', PRRS(:, :, :, 4) )
IF( LBUDGET_SV ) THEN
DO JSV = 1, NSV
CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_SV1 - 1 + JSV), 'HTURB', PRSVS(:, :, :, JSV) )
END DO
END IF
#ifdef REPRO48
#else
END IF
#endif
!----------------------------------------------------------------------------
!
!* 6. EVOLUTION OF THE TKE AND ITS DISSIPATION
! ----------------------------------------
!
! 6.1 Contribution of mass-flux in the TKE buoyancy production if
! cloud computation is not statistical
PTP = PTP + XG / PTHVREF * MZF(PFLXZTHVMF,KKA, KKU, KKL)
PTPMF=XG / PTHVREF * MZF(PFLXZTHVMF, KKA, KKU, KKL)
! 6.2 TKE evolution equation
IF (.NOT. LHARAT) THEN
!
IF (LBUDGET_TH) THEN
IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN
CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS+ ZLVOCPEXNM * PRRS(:,:,:,2) &
& + ZLSOCPEXNM * PRRS(:,:,:,4) )
ELSE IF ( KRRL >= 1 ) THEN
CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS+ ZLOCPEXNM * PRRS(:,:,:,2) )
ELSE
CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS(:, :, :) )
END IF
END IF
CALL TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKET,ZLM,ZLEPS,PDP,ZTRH, &
& PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ, &
& PTSTEP,PIMPL,ZEXPL, &
& HTURBLEN,HTURBDIM, &
& TPFILE,OTURB_DIAG, &
& PTP,PRTKES,PRTHLS,ZCOEF_DISS,PTDIFF,PTDISS,&
& TBUDGETS,KBUDGETS,&
& PEDR=PEDR)
IF (LBUDGET_TH) THEN
IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN
CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS+ ZLVOCPEXNM * PRRS(:,:,:,2) &
& + ZLSOCPEXNM * PRRS(:,:,:,4) )
ELSE IF ( KRRL >= 1 ) THEN
CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS+ ZLOCPEXNM * PRRS(:,:,:,2) )
CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS(:, :, :) )
END IF
END IF
!
ENDIF
!
!----------------------------------------------------------------------------
!
!* 7. STORES SOME INFORMATIONS RELATED TO THE TURBULENCE SCHEME
! ---------------------------------------------------------
!
IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
!
! stores the mixing length
!
TZFIELD%CMNHNAME = 'LM'
TZFIELD%CSTDNAME = ''
TZFIELD%CLONGNAME = 'LM'
TZFIELD%CUNITS = 'm'
TZFIELD%CDIR = 'XY'
TZFIELD%CCOMMENT = 'Mixing length'
TZFIELD%NGRID = 1
TZFIELD%NTYPE = TYPEREAL
TZFIELD%NDIMS = 3
TZFIELD%LTIMEDEP = .TRUE.
CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZLM)
!
IF (KRR /= 0) THEN
!
! stores the conservative potential temperature
!
TZFIELD%CMNHNAME = 'THLM'
TZFIELD%CSTDNAME = ''
TZFIELD%CLONGNAME = 'THLM'
TZFIELD%CUNITS = 'K'
TZFIELD%CDIR = 'XY'
TZFIELD%CCOMMENT = 'Conservative potential temperature'
TZFIELD%NGRID = 1
TZFIELD%NTYPE = TYPEREAL
TZFIELD%NDIMS = 3
TZFIELD%LTIMEDEP = .TRUE.
CALL IO_FIELD_WRITE(TPFILE,TZFIELD,PTHLT)
!
! stores the conservative mixing ratio
!
TZFIELD%CMNHNAME = 'RNPM'
TZFIELD%CSTDNAME = ''
TZFIELD%CLONGNAME = 'RNPM'
TZFIELD%CUNITS = 'kg kg-1'
TZFIELD%CDIR = 'XY'
TZFIELD%CCOMMENT = 'Conservative mixing ratio'
TZFIELD%NGRID = 1
TZFIELD%NTYPE = TYPEREAL
TZFIELD%NDIMS = 3
TZFIELD%LTIMEDEP = .TRUE.
CALL IO_FIELD_WRITE(TPFILE,TZFIELD,PRT(:,:,:,1))
END IF
END IF
!
!* stores value of conservative variables & wind before turbulence tendency
PDRUS_TURB = PRUS - PDRUS_TURB
PDRVS_TURB = PRVS - PDRVS_TURB
PDRTHLS_TURB = PRTHLS - PDRTHLS_TURB
PDRRTS_TURB = PRRS(:,:,:,1) - PDRRTS_TURB
PDRSVS_TURB = PRSVS - PDRSVS_TURB
!----------------------------------------------------------------------------
!
!* 8. RETRIEVE NON-CONSERVATIVE VARIABLES
! -----------------------------------
!
IF ( KRRL >= 1 ) THEN
IF ( KRRI >= 1 ) THEN
PRT(:,:,:,1) = PRT(:,:,:,1) - PRT(:,:,:,2) - PRT(:,:,:,4)
PRRS(:,:,:,1) = PRRS(:,:,:,1) - PRRS(:,:,:,2) - PRRS(:,:,:,4)
PTHLT(:,:,:) = PTHLT(:,:,:) + ZLVOCPEXNM(:,:,:) * PRT(:,:,:,2) &
+ ZLSOCPEXNM(:,:,:) * PRT(:,:,:,4)
PRTHLS(:,:,:) = PRTHLS(:,:,:) + ZLVOCPEXNM(:,:,:) * PRRS(:,:,:,2) &
+ ZLSOCPEXNM(:,:,:) * PRRS(:,:,:,4)
!
DEALLOCATE(ZLVOCPEXNM)
DEALLOCATE(ZLSOCPEXNM)
ELSE
PRT(:,:,:,1) = PRT(:,:,:,1) - PRT(:,:,:,2)
PRRS(:,:,:,1) = PRRS(:,:,:,1) - PRRS(:,:,:,2)
PTHLT(:,:,:) = PTHLT(:,:,:) + ZLOCPEXNM(:,:,:) * PRT(:,:,:,2)
PRTHLS(:,:,:) = PRTHLS(:,:,:) + ZLOCPEXNM(:,:,:) * PRRS(:,:,:,2)
END IF
END IF
! Remove non-physical negative values (unnecessary in a perfect world) + corresponding budgets
!CALL SOURCES_NEG_CORRECT(HCLOUD, 'NETUR',KRR,PTSTEP,PPABST,PTHLT,PRT,PRTHLS,PRRS,PRSVS)
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
!----------------------------------------------------------------------------
!
!* 9. LES averaged surface fluxes
! ---------------------------
!
IF (LLES_CALL) THEN
CALL SECOND_MNH(ZTIME1)
CALL LES_MEAN_SUBGRID(PSFTH,X_LES_Q0)
CALL LES_MEAN_SUBGRID(PSFRV,X_LES_E0)
DO JSV=1,NSV
CALL LES_MEAN_SUBGRID(PSFSV(:,:,JSV),X_LES_SV0(:,JSV))
END DO
CALL LES_MEAN_SUBGRID(PSFU,X_LES_UW0)
CALL LES_MEAN_SUBGRID(PSFV,X_LES_VW0)
CALL LES_MEAN_SUBGRID((PSFU*PSFU+PSFV*PSFV)**0.25,X_LES_USTAR)
!----------------------------------------------------------------------------
!
!* 10. LES for 3rd order moments
! -------------------------
!
CALL LES_MEAN_SUBGRID(ZMWTH,X_LES_SUBGRID_W2Thl)
CALL LES_MEAN_SUBGRID(ZMTH2,X_LES_SUBGRID_WThl2)
IF (KRR>0) THEN
CALL LES_MEAN_SUBGRID(ZMWR,X_LES_SUBGRID_W2Rt)
CALL LES_MEAN_SUBGRID(ZMTHR,X_LES_SUBGRID_WThlRt)
CALL LES_MEAN_SUBGRID(ZMR2,X_LES_SUBGRID_WRt2)
END IF
!
!----------------------------------------------------------------------------
!
!* 11. LES quantities depending on <w'2> in "1DIM" mode
! ------------------------------------------------
!
IF (HTURBDIM=="1DIM") THEN
CALL LES_MEAN_SUBGRID(2./3.*PTKET,X_LES_SUBGRID_U2)
X_LES_SUBGRID_V2 = X_LES_SUBGRID_U2
X_LES_SUBGRID_W2 = X_LES_SUBGRID_U2
CALL LES_MEAN_SUBGRID(2./3.*PTKET*MZF(GZ_M_W(KKA,KKU,KKL,PTHLT,PDZZ),&
KKA, KKU, KKL),X_LES_RES_ddz_Thl_SBG_W2)
IF (KRR>=1) &
CALL LES_MEAN_SUBGRID(2./3.*PTKET*MZF(GZ_M_W(KKA,KKU,KKL,PRT(:,:,:,1),PDZZ),&
&KKA, KKU, KKL),X_LES_RES_ddz_Rt_SBG_W2)
DO JSV=1,NSV
CALL LES_MEAN_SUBGRID(2./3.*PTKET*MZF(GZ_M_W(KKA,KKU,KKL,PSVT(:,:,:,JSV),PDZZ), &
&KKA, KKU, KKL), X_LES_RES_ddz_Sv_SBG_W2(:,:,:,JSV))
END DO
END IF
!----------------------------------------------------------------------------
!
!* 12. LES mixing end dissipative lengths, presso-correlations
! -------------------------------------------------------
!
CALL LES_MEAN_SUBGRID(ZLM,X_LES_SUBGRID_LMix)
CALL LES_MEAN_SUBGRID(ZLEPS,X_LES_SUBGRID_LDiss)
!
!* presso-correlations for subgrid Tke are equal to zero.
!
ZLEPS = 0. !ZLEPS is used as a work array (not used anymore)
CALL LES_MEAN_SUBGRID(ZLEPS,X_LES_SUBGRID_WP)
!
CALL SECOND_MNH(ZTIME2)
XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
END IF
!
IF(PRESENT(PLEM)) PLEM = ZLM
!----------------------------------------------------------------------------
!
IF (LHOOK) CALL DR_HOOK('TURB',1,ZHOOK_HANDLE)
CONTAINS
!
!
! ##############################################
SUBROUTINE UPDATE_ROTATE_WIND(PUSLOPE,PVSLOPE)
! ##############################################
!!
!!**** *UPDATE_ROTATE_WIND* routine to set rotate wind values at the border
!
!! AUTHOR
!! ------
!!
!! P Jabouille *CNRM METEO-FRANCE
!!
!! MODIFICATIONS
!! -------------
!! Original 24/06/99
!! J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
!!
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
USE MODE_ll
USE MODD_ARGSLIST_ll, ONLY : LIST_ll
USE MODD_CONF
!
IMPLICIT NONE
!
!* 0.1 Declarations of dummy arguments :
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: PUSLOPE,PVSLOPE
! tangential surface fluxes in the axes following the orography
!
!* 0.2 Declarations of local variables :
!
INTEGER :: IIB,IIE,IJB,IJE,IIU,IJU ! index values for the physical subdomain
TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange
INTEGER :: IINFO_ll ! return code of parallel routine
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB:UPDATE_ROTATE_WIND',0,ZHOOK_HANDLE)
!
!* 1 PROLOGUE
!
NULLIFY(TZFIELDS_ll)
CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
!
! 2 Update halo if necessary
!
!!$IF (NHALO == 1) THEN
CALL ADD2DFIELD_ll( TZFIELDS_ll, PUSLOPE, 'UPDATE_ROTATE_WIND::PUSLOPE' )
CALL ADD2DFIELD_ll( TZFIELDS_ll, PVSLOPE, 'UPDATE_ROTATE_WIND::PVSLOPE' )
CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
CALL CLEANLIST_ll(TZFIELDS_ll)
!!$ENDIF
!
! 3 Boundary conditions for non cyclic case
!
IF ( HLBCX(1) /= "CYCL" .AND. LWEST_ll()) THEN
PUSLOPE(IIB-1,:)=PUSLOPE(IIB,:)
PVSLOPE(IIB-1,:)=PVSLOPE(IIB,:)
END IF
IF ( HLBCX(2) /= "CYCL" .AND. LEAST_ll()) THEN
PUSLOPE(IIE+1,:)=PUSLOPE(IIE,:)
PVSLOPE(IIE+1,:)=PVSLOPE(IIE,:)
END IF
IF ( HLBCY(1) /= "CYCL" .AND. LSOUTH_ll()) THEN
PUSLOPE(:,IJB-1)=PUSLOPE(:,IJB)
PVSLOPE(:,IJB-1)=PVSLOPE(:,IJB)
END IF
IF( HLBCY(2) /= "CYCL" .AND. LNORTH_ll()) THEN
PUSLOPE(:,IJE+1)=PUSLOPE(:,IJE)
PVSLOPE(:,IJE+1)=PVSLOPE(:,IJE)
END IF
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
!
IF (LHOOK) CALL DR_HOOK('TURB:UPDATE_ROTATE_WIND',1,ZHOOK_HANDLE)
!
END SUBROUTINE UPDATE_ROTATE_WIND
!
! ########################################################################
SUBROUTINE COMPUTE_FUNCTION_THERMO(PALP,PBETA,PGAM,PLTT,PC,PT,PEXN,PCP,&
PLOCPEXN,PAMOIST,PATHETA )
! ########################################################################
!!
!!**** *COMPUTE_FUNCTION_THERMO* routine to compute several thermo functions
!
!! AUTHOR
!! ------
!!
!! JP Pinty *LA*
!!
!! MODIFICATIONS
!! -------------
!! Original 24/02/03
!!
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
USE MODD_CST
!
IMPLICIT NONE
!
!* 0.1 Declarations of dummy arguments
!
REAL, INTENT(IN) :: PALP,PBETA,PGAM,PLTT,PC
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
REAL, DIMENSION(:,:,:), INTENT(IN) :: PT,PEXN,PCP
!
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PLOCPEXN
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PAMOIST,PATHETA
!
!* 0.2 Declarations of local variables
!
REAL :: ZEPS ! XMV / XMD
REAL, DIMENSION(SIZE(PEXN,1),SIZE(PEXN,2),SIZE(PEXN,3)) :: ZRVSAT
REAL, DIMENSION(SIZE(PEXN,1),SIZE(PEXN,2),SIZE(PEXN,3)) :: ZDRVSATDT
!
!-------------------------------------------------------------------------------
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO',0,ZHOOK_HANDLE)
ZEPS = XMV / XMD
!
!* 1.1 Lv/Cph at t
!
PLOCPEXN(:,:,:) = ( PLTT + (XCPV-PC) * (PT(:,:,:)-XTT) ) / PCP(:,:,:)
!
!* 1.2 Saturation vapor pressure at t
!
ZRVSAT(:,:,:) = EXP( PALP - PBETA/PT(:,:,:) - PGAM*ALOG( PT(:,:,:) ) )
!
!* 1.3 saturation mixing ratio at t
!
ZRVSAT(:,:,:) = ZRVSAT(:,:,:) * ZEPS / ( PPABST(:,:,:) - ZRVSAT(:,:,:) )
!
!* 1.4 compute the saturation mixing ratio derivative (rvs')
!
ZDRVSATDT(:,:,:) = ( PBETA / PT(:,:,:) - PGAM ) / PT(:,:,:) &
* ZRVSAT(:,:,:) * ( 1. + ZRVSAT(:,:,:) / ZEPS )
!
!* 1.5 compute Amoist
!
PAMOIST(:,:,:)= 0.5 / ( 1.0 + ZDRVSATDT(:,:,:) * PLOCPEXN(:,:,:) )
!
!* 1.6 compute Atheta
!
PATHETA(:,:,:)= PAMOIST(:,:,:) * PEXN(:,:,:) * &
( ( ZRVSAT(:,:,:) - PRT(:,:,:,1) ) * PLOCPEXN(:,:,:) / &
( 1. + ZDRVSATDT(:,:,:) * PLOCPEXN(:,:,:) ) * &
( &
ZRVSAT(:,:,:) * (1. + ZRVSAT(:,:,:)/ZEPS) &
* ( -2.*PBETA/PT(:,:,:) + PGAM ) / PT(:,:,:)**2 &
+ZDRVSATDT(:,:,:) * (1. + 2. * ZRVSAT(:,:,:)/ZEPS) &
* ( PBETA/PT(:,:,:) - PGAM ) / PT(:,:,:) &
) &
- ZDRVSATDT(:,:,:) &
)
!
!* 1.7 Lv/Cph/Exner at t-1
!
PLOCPEXN(:,:,:) = PLOCPEXN(:,:,:) / PEXN(:,:,:)
!
IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO',1,ZHOOK_HANDLE)
END SUBROUTINE COMPUTE_FUNCTION_THERMO
!
! ####################
SUBROUTINE DELT(PLM,ODZ)
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
! ####################
!!
!!**** *DELT* routine to compute mixing length for DELT case
!
!! AUTHOR
!! ------
!!
!! M Tomasini *Meteo-France
!!
!! MODIFICATIONS
!! -------------
!! Original 01/05
!!
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
!* 0.1 Declarations of dummy arguments
!
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PLM
LOGICAL, INTENT(IN) :: ODZ
!
!* 0.2 Declarations of local variables
!
REAL :: ZD ! distance to the surface
!
!-------------------------------------------------------------------------------
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB:DELT',0,ZHOOK_HANDLE)
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
IF (ODZ) THEN
! Dz is take into account in the computation
DO JK = IKTB,IKTE ! 1D turbulence scheme
PLM(:,:,JK) = PZZ(:,:,JK+KKL) - PZZ(:,:,JK)
END DO
PLM(:,:,KKU) = PLM(:,:,IKE)
PLM(:,:,KKA) = PZZ(:,:,IKB) - PZZ(:,:,KKA)
IF ( HTURBDIM /= '1DIM' ) THEN ! 3D turbulence scheme
IF ( L2D) THEN
PLM(:,:,:) = SQRT( PLM(:,:,:)*MXF(PDXX(:,:,:)) )
ELSE
PLM(:,:,:) = (PLM(:,:,:)*MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./3.)
END IF
END IF
ELSE
! Dz not taken into account in computation to assure invariability with vertical grid mesh
PLM=1.E10
IF ( HTURBDIM /= '1DIM' ) THEN ! 3D turbulence scheme
IF ( L2D) THEN
PLM(:,:,:) = MXF(PDXX(:,:,:))
ELSE
PLM(:,:,:) = (MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./2.)
END IF
END IF
END IF
!
! mixing length limited by the distance normal to the surface
! (with the same factor as for BL89)
!
IF (.NOT. ORMC01) THEN
ZALPHA=0.5**(-1.5)
!
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
DO JJ=1,SIZE(PUT,2)
DO JI=1,SIZE(PUT,1)
IF (LOCEAN) THEN
DO JK=IKTE,IKTB,-1
ZD=ZALPHA*(PZZ(JI,JJ,IKTE+1)-PZZ(JI,JJ,JK))
IF ( PLM(JI,JJ,JK)>ZD) THEN
PLM(JI,JJ,JK)=ZD
ELSE
EXIT
ENDIF
END DO
ELSE
DO JK=IKTB,IKTE
ZD=ZALPHA*(0.5*(PZZ(JI,JJ,JK)+PZZ(JI,JJ,JK+KKL))&
-PZZ(JI,JJ,IKB)) *PDIRCOSZW(JI,JJ)
IF ( PLM(JI,JJ,JK)>ZD) THEN
PLM(JI,JJ,JK)=ZD
ELSE
EXIT
ENDIF
END DO
ENDIF
END DO
END DO
END IF
!
PLM(:,:,KKA) = PLM(:,:,IKB )
PLM(:,:,KKU ) = PLM(:,:,IKE)
!
IF (LHOOK) CALL DR_HOOK('TURB:DELT',1,ZHOOK_HANDLE)
END SUBROUTINE DELT
!
! ####################
SUBROUTINE DEAR(PLM)
! ####################
!!
!!**** *DEAR* routine to compute mixing length for DEARdorff case
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
!
!! AUTHOR
!! ------
!!
!! M Tomasini *Meteo-France
!!
!! MODIFICATIONS
!! -------------
!! Original 01/05
!! I.Sandu (Sept.2006) : Modification of the stability criterion
!! (theta_v -> theta_l)
!!
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
!* 0.1 Declarations of dummy arguments
!
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PLM
!
!* 0.2 Declarations of local variables
!
REAL :: ZD ! distance to the surface
REAL :: ZVAR ! Intermediary variable
REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: ZWORK2D
REAL, DIMENSION(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)) :: &
ZDTHLDZ,ZDRTDZ, &!dtheta_l/dz, drt_dz used for computing the stablity
! ! criterion
ZETHETA,ZEMOIST !coef ETHETA and EMOIST
!----------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
!
! initialize the mixing length with the mesh grid
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB:DEAR',0,ZHOOK_HANDLE)
! 1D turbulence scheme
PLM(:,:,IKTB:IKTE) = PZZ(:,:,IKTB+KKL:IKTE+KKL) - PZZ(:,:,IKTB:IKTE)
PLM(:,:,KKU) = PLM(:,:,IKE)
PLM(:,:,KKA) = PZZ(:,:,IKB) - PZZ(:,:,KKA)
IF ( HTURBDIM /= '1DIM' ) THEN ! 3D turbulence scheme
IF ( L2D) THEN
PLM(:,:,:) = SQRT( PLM(:,:,:)*MXF(PDXX(:,:,:)) )
ELSE
PLM(:,:,:) = (PLM(:,:,:)*MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./3.)
END IF
END IF
! compute a mixing length limited by the stability
!
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
ZETHETA(:,:,:) = ETHETA(KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZATHETA,PSRCT)
ZEMOIST(:,:,:) = EMOIST(KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZAMOIST,PSRCT)
!
IF (KRR>0) THEN
DO JK = IKTB+1,IKTE-1
DO JJ=1,SIZE(PUT,2)
DO JI=1,SIZE(PUT,1)
ZDTHLDZ(JI,JJ,JK)= 0.5*((PTHLT(JI,JJ,JK+KKL)-PTHLT(JI,JJ,JK ))/PDZZ(JI,JJ,JK+KKL)+ &
(PTHLT(JI,JJ,JK )-PTHLT(JI,JJ,JK-KKL))/PDZZ(JI,JJ,JK ))
ZDRTDZ(JI,JJ,JK) = 0.5*((PRT(JI,JJ,JK+KKL,1)-PRT(JI,JJ,JK ,1))/PDZZ(JI,JJ,JK+KKL)+ &
(PRT(JI,JJ,JK ,1)-PRT(JI,JJ,JK-KKL,1))/PDZZ(JI,JJ,JK ))
IF (LOCEAN) THEN
ZVAR=XG*(XALPHAOC*ZDTHLDZ(JI,JJ,JK)-XBETAOC*ZDRTDZ(JI,JJ,JK))
ELSE
ZVAR=XG/PTHVREF(JI,JJ,JK)* &
(ZETHETA(JI,JJ,JK)*ZDTHLDZ(JI,JJ,JK)+ZEMOIST(JI,JJ,JK)*ZDRTDZ(JI,JJ,JK))
END IF
!
IF (ZVAR>0.) THEN
PLM(JI,JJ,JK)=MAX(XMNH_EPSILON,MIN(PLM(JI,JJ,JK), &
0.76* SQRT(PTKET(JI,JJ,JK)/ZVAR)))
END IF
END DO
END DO
END DO
ELSE! For dry atmos or unsalted ocean runs
DO JK = IKTB+1,IKTE-1
DO JJ=1,SIZE(PUT,2)
DO JI=1,SIZE(PUT,1)
ZDTHLDZ(JI,JJ,JK)= 0.5*((PTHLT(JI,JJ,JK+KKL)-PTHLT(JI,JJ,JK ))/PDZZ(JI,JJ,JK+KKL)+ &
(PTHLT(JI,JJ,JK )-PTHLT(JI,JJ,JK-KKL))/PDZZ(JI,JJ,JK ))
IF (LOCEAN) THEN
ZVAR= XG*XALPHAOC*ZDTHLDZ(JI,JJ,JK)
ELSE
ZVAR= XG/PTHVREF(JI,JJ,JK)*ZETHETA(JI,JJ,JK)*ZDTHLDZ(JI,JJ,JK)
END IF
IF (ZVAR>0.) THEN
PLM(JI,JJ,JK)=MAX(XMNH_EPSILON,MIN(PLM(JI,JJ,JK), &
0.76* SQRT(PTKET(JI,JJ,JK)/ZVAR)))
END IF
END DO
END DO
END DO
END IF
! special case near the surface
ZDTHLDZ(:,:,IKB)=(PTHLT(:,:,IKB+KKL)-PTHLT(:,:,IKB))/PDZZ(:,:,IKB+KKL)
! For dry simulations
IF (KRR>0) THEN
ZDRTDZ(:,:,IKB)=(PRT(:,:,IKB+KKL,1)-PRT(:,:,IKB,1))/PDZZ(:,:,IKB+KKL)
ELSE
ZDRTDZ(:,:,IKB)=0
ENDIF
IF (LOCEAN) THEN
ZWORK2D(:,:)=XG*(XALPHAOC*ZDTHLDZ(:,:,IKB)-XBETAOC*ZDRTDZ(:,:,IKB))
ELSE
ZWORK2D(:,:)=XG/PTHVREF(:,:,IKB)* &
(ZETHETA(:,:,IKB)*ZDTHLDZ(:,:,IKB)+ZEMOIST(:,:,IKB)*ZDRTDZ(:,:,IKB))
END IF
WHERE(ZWORK2D(:,:)>0.)
PLM(:,:,IKB)=MAX(XMNH_EPSILON,MIN( PLM(:,:,IKB), &
0.76* SQRT(PTKET(:,:,IKB)/ZWORK2D(:,:))))
END WHERE
!
! mixing length limited by the distance normal to the surface (with the same factor as for BL89)
!
IF (.NOT. ORMC01) THEN
ZALPHA=0.5**(-1.5)
!
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
DO JJ=1,SIZE(PUT,2)
DO JI=1,SIZE(PUT,1)
IF (LOCEAN) THEN
DO JK=IKTE,IKTB,-1
ZD=ZALPHA*(PZZ(JI,JJ,IKTE+1)-PZZ(JI,JJ,JK))
IF ( PLM(JI,JJ,JK)>ZD) THEN
PLM(JI,JJ,JK)=ZD
ELSE
EXIT
ENDIF
END DO
ELSE
DO JK=IKTB,IKTE
ZD=ZALPHA*(0.5*(PZZ(JI,JJ,JK)+PZZ(JI,JJ,JK+KKL))-PZZ(JI,JJ,IKB)) &
*PDIRCOSZW(JI,JJ)
IF ( PLM(JI,JJ,JK)>ZD) THEN
PLM(JI,JJ,JK)=ZD
ELSE
EXIT
ENDIF
END DO
ENDIF
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
END DO
END DO
END IF
!
PLM(:,:,KKA) = PLM(:,:,IKB )
PLM(:,:,IKE ) = PLM(:,:,IKE-KKL)
PLM(:,:,KKU ) = PLM(:,:,KKU-KKL)
!
IF (LHOOK) CALL DR_HOOK('TURB:DEAR',1,ZHOOK_HANDLE)
END SUBROUTINE DEAR
!
! #########################
SUBROUTINE CLOUD_MODIF_LM
! #########################
!!
!!*****CLOUD_MODIF_LM routine to:
!! 1/ change the mixing length in the clouds
!! 2/ emphasize the mixing length in the cloud
!! by the coefficient ZCOEF_AMPL calculated here
!! when the CEI index is above ZCEI_MIN.
!!
!!
!! ZCOEF_AMPL ^
!! |
!! |
!! ZCOEF_AMPL_SAT - ---------- Saturation
!! (XDUMMY1) | -
!! | -
!! | -
!! | -
!! | - Amplification
!! | - straight
!! | - line
!! | -
!! | -
!! | -
!! | -
!! | -
!! 1 ------------
!! |
!! |
!! 0 -----------|------------|----------> PCEI
!! 0 ZCEI_MIN ZCEI_MAX
!! (XDUMMY2) (XDUMMY3)
!!
!!
!!
!! AUTHOR
!! ------
!! M. Tomasini *CNRM METEO-FRANCE
!!
!! MODIFICATIONS
!! -------------
!! Original 09/07/04
!!
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
IMPLICIT NONE
!
REAL :: ZPENTE ! Slope of the amplification straight line
REAL :: ZCOEF_AMPL_CEI_NUL! Ordonnate at the origin of the
! amplification straight line
REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZCOEF_AMPL
! Amplification coefficient of the mixing length
! when the instability criterium is verified
REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZLM_CLOUD
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
! Turbulent mixing length in the clouds
!
!-------------------------------------------------------------------------------
!
!* 1. INITIALISATION
! --------------
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB:CLOUD_MODIF_LM',0,ZHOOK_HANDLE)
ZPENTE = ( PCOEF_AMPL_SAT - 1. ) / ( PCEI_MAX - PCEI_MIN )
ZCOEF_AMPL_CEI_NUL = 1. - ZPENTE * PCEI_MIN
!
ZCOEF_AMPL(:,:,:) = 1.
!
!* 2. CALCULATION OF THE AMPLIFICATION COEFFICIENT
! --------------------------------------------
!
! Saturation
!
WHERE ( PCEI(:,:,:)>=PCEI_MAX ) ZCOEF_AMPL(:,:,:)=PCOEF_AMPL_SAT
!
! Between the min and max limits of CEI index, linear variation of the
! amplification coefficient ZCOEF_AMPL as a function of CEI
!
WHERE ( PCEI(:,:,:) < PCEI_MAX .AND. &
PCEI(:,:,:) > PCEI_MIN ) &
ZCOEF_AMPL(:,:,:) = ZPENTE * PCEI(:,:,:) + ZCOEF_AMPL_CEI_NUL
!
!
!* 3. CALCULATION OF THE MIXING LENGTH IN CLOUDS
! ------------------------------------------
!
IF (HTURBLEN_CL == HTURBLEN) THEN
ZLM_CLOUD(:,:,:) = ZLM(:,:,:)
ELSE
SELECT CASE (HTURBLEN_CL)
!
!* 3.1 BL89 mixing length
! ------------------
CASE ('BL89','RM17','ADAP')
CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM_CLOUD)
!
!* 3.2 Delta mixing length
! -------------------
CASE ('DELT')
CALL DELT(ZLM_CLOUD,ODZ=.TRUE.)
!
!* 3.3 Deardorff mixing length
! -----------------------
CASE ('DEAR')
CALL DEAR(ZLM_CLOUD)
!
END SELECT
ENDIF
!
!* 4. MODIFICATION OF THE MIXING LENGTH IN THE CLOUDS
! -----------------------------------------------
!
! Impression before modification of the mixing length
IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
TZFIELD%CMNHNAME = 'LM_CLEAR_SKY'
TZFIELD%CSTDNAME = ''
TZFIELD%CLONGNAME = 'LM_CLEAR_SKY'
TZFIELD%CUNITS = 'm'
TZFIELD%CDIR = 'XY'
TZFIELD%CCOMMENT = 'X_Y_Z_LM CLEAR SKY'
TZFIELD%NGRID = 1
TZFIELD%NTYPE = TYPEREAL
TZFIELD%NDIMS = 3
TZFIELD%LTIMEDEP = .TRUE.
CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZLM)
ENDIF
!
! Amplification of the mixing length when the criteria are verified
!
WHERE (ZCOEF_AMPL(:,:,:) /= 1.) ZLM(:,:,:) = ZCOEF_AMPL(:,:,:)*ZLM_CLOUD(:,:,:)
!
! Cloud mixing length in the clouds at the points which do not verified the CEI
!
WHERE (PCEI(:,:,:) == -1.) ZLM(:,:,:) = ZLM_CLOUD(:,:,:)
!
!
!* 5. IMPRESSION
! ----------
!
IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
TZFIELD%CMNHNAME = 'COEF_AMPL'
TZFIELD%CSTDNAME = ''
TZFIELD%CLONGNAME = 'COEF_AMPL'
TZFIELD%CUNITS = '1'
TZFIELD%CDIR = 'XY'
TZFIELD%CCOMMENT = 'X_Y_Z_COEF AMPL'
TZFIELD%NGRID = 1
TZFIELD%NTYPE = TYPEREAL
TZFIELD%NDIMS = 3
TZFIELD%LTIMEDEP = .TRUE.
CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZCOEF_AMPL)
TZFIELD%CMNHNAME = 'LM_CLOUD'
TZFIELD%CSTDNAME = ''
TZFIELD%CLONGNAME = 'LM_CLOUD'
TZFIELD%CUNITS = 'm'
TZFIELD%CDIR = 'XY'
TZFIELD%CCOMMENT = 'X_Y_Z_LM CLOUD'
TZFIELD%NGRID = 1
TZFIELD%NTYPE = TYPEREAL
TZFIELD%NDIMS = 3
CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZLM_CLOUD)
!
ENDIF
!
IF (LHOOK) CALL DR_HOOK('TURB:CLOUD_MODIF_LM',1,ZHOOK_HANDLE)
END SUBROUTINE CLOUD_MODIF_LM
!
END SUBROUTINE TURB