diff --git a/src/MNH/diag.f90 b/src/MNH/diag.f90 index 7324f3dd187bb7e1def5a8358913ee998d6a87d7..af74d5cec8ea88088f53880134d179acd3df0f56 100644 --- a/src/MNH/diag.f90 +++ b/src/MNH/diag.f90 @@ -81,6 +81,7 @@ !! 10/2017 (G.Delautier) New boundary layer height : replace LBLTOP by CBLTOP !! 10/2017 (T Dauhut) Add parallel 3D clustering !! 01/2018 (J.-P. Chaboureau) Add altitude interpolation +!! 01/2018 (J.-P. Chaboureau) Add coarse graining !------------------------------------------------------------------------------- ! !* 0. DECLARATIONS @@ -230,7 +231,8 @@ NAMELIST/NAM_DIAG/ CISO, LVAR_RS, LVAR_LS, & XGRID,NBELEV,XELEV,NBRAD,LQUAD,LFALL,LWBSCS,LWREFL,& XREFLMIN,XREFLVDOPMIN,LSNRT,XSNRMIN,& LLIDAR,CVIEW_LIDAR,XALT_LIDAR,XWVL_LIDAR,& - LISOPR,XISOPR,LISOTH,XISOTH,LISOAL,XISOAL,LHU_FLX,LVISI,LLIMA_DIAG,& + LISOPR,XISOPR,LISOTH,XISOTH,LISOAL,XISOAL,LCOARSE,NDXCOARSE, & + LHU_FLX,LVISI,LLIMA_DIAG,& LCLSTR,LBOTUP,CFIELD,XTHRES ! NAMELIST/NAM_DIAG_FILE/ YINIFILE,YINIFILEPGD, YSUFFIX @@ -384,6 +386,9 @@ XISOTH(:)=0. LISOAL=.FALSE. XISOAL(:)=-1. ! +LCOARSE=.FALSE. +NDXCOARSE=1 +! !------------------------------------------------------------------------------- ! !* 1.0 Namelist reading diff --git a/src/MNH/modd_diag_flag.f90 b/src/MNH/modd_diag_flag.f90 index 118897cadf70197ed6b61f90f9a5f28ffada1d4c..f5ac8fbcf552809f68b528fab0fc8ff86ef78b6c 100644 --- a/src/MNH/modd_diag_flag.f90 +++ b/src/MNH/modd_diag_flag.f90 @@ -42,6 +42,7 @@ !! 10/2017 (G.Delautier) New boundary layer height : replace LBLTOP by CBLTOP !! T. Dauhut 10/2017 add parallel 3D clustering !! J.-P. Chaboureau 01/2018 add altitude interpolation +!! J.-P. Chaboureau 01/2018 add coarse graining !! !------------------------------------------------------------------------------- ! @@ -138,6 +139,9 @@ REAL,DIMENSION(10):: XISOTH ! list of level for isentropic interpolation LOGICAL :: LISOAL ! flag to write on altitude level REAL,DIMENSION(99):: XISOAL ! list of level for altitude interpolation ! +LOGICAL :: LCOARSE ! flag for coarse graining +INTEGER :: NDXCOARSE ! gridpoint number for coarse graining +! LOGICAL :: LCLSTR ! flag for 3D clustering LOGICAL :: LBOTUP ! to propagate clustering from bottom to top, otherwise top to bottom CHARACTER(LEN=8) :: CFIELD ! field on which clustering is applied, could be 'W' or 'CLOUD' diff --git a/src/MNH/neighboravg.f90 b/src/MNH/neighboravg.f90 new file mode 100644 index 0000000000000000000000000000000000000000..903d0ad280a81063f41477bc28e959b2a46f91d4 --- /dev/null +++ b/src/MNH/neighboravg.f90 @@ -0,0 +1,253 @@ +!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier +!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence +!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!MNH_LIC for details. version 1. +!###################### +MODULE MODE_NEIGHBORAVG +!###################### +! +CONTAINS +! +SUBROUTINE BLOCKAVG(PMATIN,KDX,KDY,PMATOUT) +! ##################################################################### +! +!! PURPOSE +!! ------- +!! AVERAGE PMATIN FIELDS BY BLOCK OF DX TIMES DY +!! +!! AUTHOR +!! ------ +!! J. Escobar *L.A.* +!! +!! MODIFICATIONS +!! ------------- +!! +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +USE MODD_PARAMETERS, ONLY : JPHEXT +USE MODD_LBC_N , ONLY : CLBCX, CLBCY +USE MODE_ll , ONLY : GET_INDICE_ll, GET_OR_ll +USE MODE_ll , ONLY : LNORTH_ll, LSOUTH_ll, LEAST_ll, LWEST_ll +USE MODI_GET_HALO , ONLY : GET_HALO + +IMPLICIT NONE +! +REAL, INTENT(IN),DIMENSION(:,:,:) :: PMATIN +INTEGER, INTENT(IN) :: KDX,KDY +REAL, INTENT(INOUT),DIMENSION(:,:,:) :: PMATOUT +! +!* 0.1 declarations of local variables +! +REAL, DIMENSION(SIZE(PMATIN,1),SIZE(PMATIN,2),SIZE(PMATIN,3)) :: ZSUM, ZTMP + +INTEGER :: IIB, IIE, IJB, IJE +INTEGER :: IIBDX, IJBDY, IDX, IDY +INTEGER :: IDIM1, IDIM2, IDIM3 +INTEGER :: II, IJ, IXOR_ll, IYOR_ll +!------------------------------------------------------------------------------- +! +!* 1. INITIALIZATION +! -------------- +! +! Init the local grid info + +CALL GET_INDICE_ll (IIB,IJB,IIE,IJE) +CALL GET_OR_ll('B', IXOR_ll, IYOR_ll) + +! shift of the first local pts to DX/DY grid pts + +IDX = MOD(IIB+IXOR_ll-1-JPHEXT-1,KDX) +IDY = MOD(IJB+IYOR_ll-1-JPHEXT-1,KDY) + +IF ( IDX .NE. 0 ) IDX = KDX - IDX +IF ( IDY .NE. 0 ) IDY = KDY - IDY + +IIBDX = IIB + IDX +IJBDY = IJB + IDY + +! Init Sum + +CALL GET_HALO(PMATIN) +ZSUM = 0.0 +ZTMP = PMATIN + +ZSUM(IIBDX:IIE:KDX,:,:) = PMATIN(IIBDX:IIE:KDX,:,:) + +! Do the sum on mod(X,KDX) column + +DO II = 1, KDX-1 + IF (LEAST_ll() .AND. CLBCX(2)/='CYCL') ZTMP(IIE+1,:,:) = 0.0 + CALL GET_HALO(ZTMP) + ZTMP(IIB:IIE,:,:) = ZTMP(IIB+1:IIE+1,:,:) + ZSUM(IIBDX:IIE:KDX,:,:) = ZSUM(IIBDX:IIE:KDX,:,:) + ZTMP(IIBDX:IIE:KDX,:,:) +END DO + +ZTMP = ZSUM + +! DO the sum on mod(Y,KDY) raw + +DO IJ = 1, KDY-1 + IF (LNORTH_ll() .AND. CLBCY(2)/='CYCL') ZTMP(:,IJE+1,:) = 0.0 + CALL GET_HALO(ZTMP) + ZTMP(IIB:IIE,IJB:IJE,:) = ZTMP(IIB:IIE,IJB+1:IJE+1,:) + ZSUM(IIBDX:IIE:KDX,IJBDY:IJE:KDY,:) = ZSUM(IIBDX:IIE:KDX,IJBDY:IJE:KDY,:) & + + ZTMP(IIBDX:IIE:KDX,IJBDY:IJE:KDY,:) +END DO + +ZTMP = 0.0 + +! Dispatch sum on mod(Y,KDY) raw + +ZTMP(IIBDX:IIE:KDX,IJBDY:IJE:KDY,:) = ZSUM(IIBDX:IIE:KDX,IJBDY:IJE:KDY,:) + +DO IJ = 1, KDY-1 + IF (LSOUTH_ll() .AND. CLBCY(1)/='CYCL') ZTMP(:,IJB-1,:) = 0.0 + CALL GET_HALO(ZTMP) + ZTMP(IIB:IIE,IJB:IJE,:) = ZTMP(IIB:IIE,IJB-1:IJE-1,:) + ZTMP(IIBDX:IIE:KDX,IJBDY:IJE:KDY,:) = ZSUM(IIBDX:IIE:KDX,IJBDY:IJE:KDY,:) +END DO + + +! Dispatch sum on mod(X,KDX) column + +ZSUM(IIBDX:IIE:KDX,:,:) = ZTMP(IIBDX:IIE:KDX,:,:) + +DO II = 1, KDX-1 + IF (LWEST_ll() .AND. CLBCX(1)/='CYCL') ZTMP(IIB-1,:,:) = 0.0 + CALL GET_HALO(ZTMP) + ZTMP(IIB:IIE,:,:) = ZTMP(IIB-1:IIE-1,:,:) + ZTMP(IIBDX:IIE:KDX,:,:) = ZSUM(IIBDX:IIE:KDX,:,:) +END DO + +CALL GET_HALO(ZTMP) +PMATOUT(:,:,:) = ZTMP(:,:,:) / float(KDX*KDY) + +END SUBROUTINE BLOCKAVG + +SUBROUTINE MOVINGAVG(PMATIN,KDX,KDY,PMATOUT) +! ##################################################################### +! +!! PURPOSE +!! ------- +!! MOVING AVERAGE PMATIN FIELDS OVER OF DX TIMES DY +!! +!! AUTHOR +!! ------ +!! J. Escobar *L.A.* +!! +!! MODIFICATIONS +!! ------------- +!! +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +! +USE MODD_PARAMETERS, ONLY : JPHEXT +USE MODD_LBC_N , ONLY : CLBCX,CLBCY +USE MODE_ll , ONLY : GET_INDICE_ll , GET_OR_ll +USE MODE_ll , ONLY : LNORTH_ll , LSOUTH_ll, LEAST_ll , LWEST_ll +USE MODI_GET_HALO , ONLY : GET_HALO + +IMPLICIT NONE +! +REAL, INTENT(IN),DIMENSION(:,:,:) :: PMATIN +INTEGER, INTENT(IN) :: KDX,KDY +REAL, INTENT(INOUT),DIMENSION(:,:,:) :: PMATOUT +!local var + +REAL, DIMENSION(SIZE(PMATIN,1),SIZE(PMATIN,2),SIZE(PMATIN,3)) :: ZSUMP1 , ZSUMM1 , ZTMP , ZTMP2 + +INTEGER :: IIB,IIE,IJB,IJE + +INTEGER :: II, IJ, IXOR_ll, IYOR_ll , ISX + + ! Init the local grid info + +CALL GET_INDICE_ll (IIB,IJB,IIE,IJE) + + ! Init Sum + +ZTMP = PMATIN + + ! shift input tab to -KDX , -KDY + +CALL SHIFTXY(ZTMP,-KDX,-KDY) + +ZSUMP1 = 0.0 + +ISX = +1 +DO IJ = 1 , 2*KDY +1 +! Do the sum on X+/-1 + DO II = 1 , 2*KDX +1 + ZSUMP1 = ZSUMP1 + ZTMP + IF ( II .NE. 2*KDY +1 ) THEN + CALL SHIFTXY(ZTMP,ISX,0) + END IF + END DO +! Do the sum on Y+1 + CALL SHIFTXY(ZTMP,0,1) + ISX = - ISX +END DO + +PMATOUT(:,:,:) = ZSUMP1(:,:,:) / FLOAT((1+2*KDX)*(1+2*KDY)) + +END SUBROUTINE MOVINGAVG + +SUBROUTINE SHIFTXY(PTAB,KDX,KDY) + +USE MODE_ll , ONLY : GET_INDICE_ll , GET_OR_ll +USE MODE_ll , ONLY : LNORTH_ll , LSOUTH_ll, LEAST_ll , LWEST_ll +USE MODI_GET_HALO , ONLY : GET_HALO + +USE MODD_LBC_N , ONLY : CLBCX,CLBCY + +IMPLICIT NONE + +!----- Argument + +REAL ,DIMENSION(:,:,:) :: PTAB +INTEGER, INTENT(IN) :: KDX,KDY + +!----- Local var + +INTEGER :: IIB,IIE,IJB,IJE +INTEGER :: ISDX,ISDY,IADX,IADY +INTEGER :: II, IJ + +!------ init + +CALL GET_INDICE_ll (IIB,IJB,IIE,IJE) + +ISDX = SIGN(1,KDX) +ISDY = SIGN(1,KDY) +IADX = ABS(KDX) +IADY = ABS(KDY) + +!---------- shift in X +DO II = 1 , IADX + IF (LEAST_ll() .AND. CLBCX(2)/='CYCL' ) PTAB(IIE+1,:,:) = 0.0 + IF (LWEST_ll() .AND. CLBCX(1)/='CYCL' ) PTAB(IIB-1,:,:) = 0.0 + CALL GET_HALO(PTAB) + PTAB(IIB:IIE,:,:) = PTAB(IIB+ISDX:IIE+ISDX,:,:) +END DO + +!---------- shitf in Y + +DO IJ = 1 , IADY + IF (LNORTH_ll() .AND. CLBCY(2)/='CYCL') PTAB(:,IJE+1,:) = 0.0 + IF (LSOUTH_ll() .AND. CLBCY(1)/='CYCL') PTAB(:,IJE+1,:) = 0.0 + CALL GET_HALO(PTAB) + PTAB(IIB:IIE,IJB:IJE,:) = PTAB(IIB:IIE,IJB+ISDY:IJE+ISDY,:) +END DO + +CALL GET_HALO(PTAB) + +END SUBROUTINE SHIFTXY + +END MODULE MODE_NEIGHBORAVG + diff --git a/src/MNH/write_lfifm1_for_diag_supp.f90 b/src/MNH/write_lfifm1_for_diag_supp.f90 index b4b8291860b70dce21a07e29d9025e31c3805ec3..2a778a69bbfefe3a25189171770efd3bfae31db7 100644 --- a/src/MNH/write_lfifm1_for_diag_supp.f90 +++ b/src/MNH/write_lfifm1_for_diag_supp.f90 @@ -87,6 +87,7 @@ END MODULE MODI_WRITE_LFIFM1_FOR_DIAG_SUPP !! F. Brosse 10/2016 add chemical production destruction terms outputs !! M.Leriche 01/07/2017 Add DIAG chimical surface fluxes !! J.-P. Chaboureau 01/2018 add altitude interpolation +!! J.-P. Chaboureau 01/2018 add coarse graining !------------------------------------------------------------------------------- ! !* 0. DECLARATIONS @@ -142,6 +143,7 @@ USE MODI_GRADIENT_V USE MODI_GRADIENT_UV ! USE MODI_SHUMAN +USE MODE_NEIGHBORAVG #ifdef MNH_RTTOV_8 USE MODI_CALL_RTTOV8 #endif @@ -216,6 +218,12 @@ INTEGER :: IAL REAL :: ZFILLVAL REAL, DIMENSION(:), ALLOCATABLE :: ZAL REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZWAL +! +! variables needed for coarse graining +REAL,DIMENSION(SIZE(XTHT,1),SIZE(XTHT,2),SIZE(XTHT,3)) :: ZUT_PRM,ZVT_PRM,ZWT_PRM +REAL,DIMENSION(SIZE(XTHT,1),SIZE(XTHT,2),SIZE(XTHT,3)) :: ZUU_AVG,ZVV_AVG,ZWW_AVG +INTEGER :: IDX +CHARACTER(LEN=3) :: YDX !------------------------------------------------------------------------------- ! !* 0. ARRAYS BOUNDS INITIALIZATION @@ -1250,7 +1258,63 @@ END IF ! !------------------------------------------------------------------------------- ! -!* 7. DIAGNOSTIC RELATED TO CHEMISTRY +!* 7. COARSE GRAINING DIAGNOSTIC +! -------------------------- +! +IF (LCOARSE) THEN + IDX = NDXCOARSE +!------------------------------- +! AVERAGE OF TKE BY BLOCK OF IDX POINTS + CALL BLOCKAVG(XUT,IDX,IDX,ZWORK31) + ZUT_PRM=XUT-ZWORK31 + CALL BLOCKAVG(XVT,IDX,IDX,ZWORK31) + ZVT_PRM=XVT-ZWORK31 + CALL BLOCKAVG(XWT,IDX,IDX,ZWORK31) + ZWT_PRM=XWT-ZWORK31 +! + ZWORK31=MXF(ZUT_PRM*ZUT_PRM) + CALL BLOCKAVG(ZWORK31,IDX,IDX,ZUU_AVG) + ZWORK31=MYF(ZVT_PRM*ZVT_PRM) + CALL BLOCKAVG(ZWORK31,IDX,IDX,ZVV_AVG) + ZWORK31=MZF(1,IKU,1,ZWT_PRM*ZWT_PRM) + CALL BLOCKAVG(ZWORK31,IDX,IDX,ZWW_AVG) + CALL BLOCKAVG(XTKET,IDX,IDX,ZWORK31) + ZWORK31=0.5*( ZUU_AVG+ZVV_AVG+ZWW_AVG ) + ZWORK31 + WRITE (YDX,FMT='(I3.3)') IDX + YRECFM='TKEBAVG'//YDX + YCOMMENT='TKE_BLOCKAVG'//YDX + IGRID = 1 + ILENCH = LEN(YCOMMENT) + CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,'XY',ZWORK31,IGRID,ILENCH,YCOMMENT,IRESP) +!--------------------------------- +! MOVING AVERAGE OF TKE OVER IDX+1 POINTS + IDX = IDX/2 + CALL MOVINGAVG(XUT,IDX,IDX,ZWORK31) + ZUT_PRM=XUT-ZWORK31 + CALL MOVINGAVG(XVT,IDX,IDX,ZWORK31) + ZVT_PRM=XVT-ZWORK31 + CALL MOVINGAVG(XWT,IDX,IDX,ZWORK31) + ZWT_PRM=XWT-ZWORK31 +! + ZWORK31=MXF(ZUT_PRM*ZUT_PRM) + CALL MOVINGAVG(ZWORK31,IDX,IDX,ZUU_AVG) + ZWORK31=MYF(ZVT_PRM*ZVT_PRM) + CALL MOVINGAVG(ZWORK31,IDX,IDX,ZVV_AVG) + ZWORK31=MZF(1,IKU,1,ZWT_PRM*ZWT_PRM) + CALL MOVINGAVG(ZWORK31,IDX,IDX,ZWW_AVG) + CALL MOVINGAVG(XTKET,IDX,IDX,ZWORK31) + ZWORK31=0.5*( ZUU_AVG+ZVV_AVG+ZWW_AVG ) + ZWORK31 + WRITE (YDX,FMT='(I3.3)') 2*IDX+1 + YRECFM='TKEMAVG'//YDX + YCOMMENT='TKE_MOVINGAVG'//YDX + IGRID = 1 + ILENCH = LEN(YCOMMENT) + CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,'XY',ZWORK31,IGRID,ILENCH,YCOMMENT,IRESP) +END IF +! +!------------------------------------------------------------------------------- +! +!* 8. DIAGNOSTIC RELATED TO CHEMISTRY ! ------------------------------- ! IF (NEQ_BUDGET>0) THEN