Newer
Older
RODIER Quentin
committed
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(TZFILE,TZFIELD,ZLM)
!
IF (KRR /= 0) THEN
!
! stores the conservative potential temperature
!
RODIER Quentin
committed
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(TZFILE,TZFIELD,PTHLT)
!
! stores the conservative mixing ratio
!
RODIER Quentin
committed
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(TZFILE,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
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
!----------------------------------------------------------------------------
!
!* 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)
CALL LES_MEAN_SUBGRID(2./3.*PTKET,X_LES_SUBGRID_V2)
CALL LES_MEAN_SUBGRID(2./3.*PTKET,X_LES_SUBGRID_W2)
CALL LES_MEAN_SUBGRID(2./3.*PTKET*MZF(GZ_M_W(PTHLT,PDZZ, KKA, KKU, KKL),&
KKA, KKU, KKL),X_LES_RES_ddz_Thl_SBG_W2)
IF (KRR>=1) &
CALL LES_MEAN_SUBGRID(2./3.*PTKET*MZF(GZ_M_W(PRT(:,:,:,1),PDZZ, KKA, KKU, KKL),&
&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(PSVT(:,:,:,JSV),PDZZ, KKA, KKU, KKL), &
&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
RODIER Quentin
committed
!! J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
!!
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
RODIER Quentin
committed
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 :
!
RODIER Quentin
committed
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
!
RODIER Quentin
committed
NULLIFY(TZFIELDS_ll)
RODIER Quentin
committed
IIU=SIZE(PUSLOPE,1)
IJU=SIZE(PUSLOPE,2)
CALL GET_INDICE_ll (IIB,IJB,IIE,IJE,IIU,IJU)
!
! 2 Update halo if necessary
!
RODIER Quentin
committed
!!$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
!
RODIER Quentin
committed
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
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
!
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
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
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
!
! ####################
RODIER Quentin
committed
SUBROUTINE DELT(PLM,ODZ)
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
! ####################
!!
!!**** *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
RODIER Quentin
committed
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)
RODIER Quentin
committed
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
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)
!
DO JJ=1,SIZE(PUT,2)
DO JI=1,SIZE(PUT,1)
RODIER Quentin
committed
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)
! ####################
!!
RODIER Quentin
committed
!!**** *DEAR* routine to compute mixing length for DEARdorff case
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
!
!! 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
RODIER Quentin
committed
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)
RODIER Quentin
committed
! 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
!
ZETHETA(:,:,:) = ETHETA(KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZATHETA,PSRCT)
ZEMOIST(:,:,:) = EMOIST(KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZAMOIST,PSRCT)
RODIER Quentin
committed
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
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)
RODIER Quentin
committed
! For dry simulations
IF (KRR>0) THEN
ZDRTDZ(:,:,IKB)=(PRT(:,:,IKB+KKL,1)-PRT(:,:,IKB,1))/PDZZ(:,:,IKB+KKL)
ELSE
ZDRTDZ(:,:,IKB)=0
ENDIF
RODIER Quentin
committed
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.)
RODIER Quentin
committed
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)
!
DO JJ=1,SIZE(PUT,2)
DO JI=1,SIZE(PUT,1)
RODIER Quentin
committed
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
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
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
1625
1626
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
! 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
! ------------------
RODIER Quentin
committed
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')
RODIER Quentin
committed
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. TZFILE%LOPENED ) THEN
RODIER Quentin
committed
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(TZFILE,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. TZFILE%LOPENED ) THEN
RODIER Quentin
committed
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(TZFILE,TZFIELD,ZCOEF_AMPL)
RODIER Quentin
committed
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(TZFILE,TZFIELD,ZLM_CLOUD)
!
ENDIF
!
IF (LHOOK) CALL DR_HOOK('TURB:CLOUD_MODIF_LM',1,ZHOOK_HANDLE)
END SUBROUTINE CLOUD_MODIF_LM
!
END SUBROUTINE TURB