!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. !----------------------------------------------------------------- !--------------- special set of characters for RCS information !----------------------------------------------------------------- ! $Source: /srv/cvsroot/MNH-VX-Y-Z/src/MNH/num_diff.f90,v $ $Revision: 1.2.2.1.2.4.16.1.2.2.2.1 $ $Date: 2015/09/16 14:31:20 $ !----------------------------------------------------------------- !----------------------------------------------------------------- ! #################### MODULE MODI_NUM_DIFF2 ! #################### ! INTERFACE ! SUBROUTINE NUM_DIFF2 ( HLBCX, HLBCY, KRR, KSV,PDK2U, PDK4U, & PDK2TH, PDK4TH,PDK2SV, PDK4SV, KMI, & PUT, PVT, PWT, PTHT, PTKET, PRT, PSVT, & PLSUM,PLSVM,PLSWM,PLSTHM,PLSRVM,PRHODJ, & PRUS, PRVS, PRWS, PRTHS, PRTKES, PRRS, PRSVS, & OZDIFFU,ONUMDIFU, ONUMDIFTH, ONUMDIFSV, & TPHALO2LIST, TPHALO2LSLIST,PZDIFFU_HALO2 ,KTCOUNT ) ! USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll USE MODE_TYPE_ZDIFFU ! CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! X direction LBC type CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! Y direction LBC type ! INTEGER, INTENT(IN) :: KRR ! Number of moist variables INTEGER, INTENT(IN) :: KSV ! Number of Scalar Variables INTEGER, INTENT(IN) :: KMI ! Model index INTEGER, INTENT(IN) :: KTCOUNT ! REAL, INTENT(IN) ::PDK2U ! 2nd order dif. coef. /dx2 REAL, INTENT(IN) :: PDK4U ! 4th order dif. coef. /dx4 ! for momentum REAL, INTENT(IN) ::PDK2TH ! for theta and mixing ratios REAL, INTENT(IN) :: PDK4TH ! and TKE REAL, INTENT(IN) ::PDK2SV ! for theta and mixing ratios REAL, INTENT(IN) :: PDK4SV ! and TKE ! REAL, DIMENSION(:,:,:), INTENT(IN) :: PLSUM,PLSVM ! Larger REAL, DIMENSION(:,:,:), INTENT(IN) :: PLSWM,PLSTHM ! Scales REAL, DIMENSION(:,:,:), INTENT(IN) :: PLSRVM ! Variables ! REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry density of ! ! anelastic reference state * Jacobian ! REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT, PWT , PTKET REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT , PSVT ! REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS, PRWS ! Source terms REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTHS, PRTKES REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS , PRSVS ! ! LOGICAL, INTENT(IN) :: OZDIFFU LOGICAL, INTENT(IN) :: ONUMDIFU LOGICAL, INTENT(IN) :: ONUMDIFTH LOGICAL, INTENT(IN) :: ONUMDIFSV TYPE(HALO2LIST_ll), POINTER :: TPHALO2LIST, TPHALO2LSLIST ! list for diffusion TYPE(TYPE_ZDIFFU_HALO2) :: PZDIFFU_HALO2 ! END SUBROUTINE NUM_DIFF2 ! END INTERFACE ! END MODULE MODI_NUM_DIFF2 ! ! ! ! ###################################################################### SUBROUTINE NUM_DIFF2 ( HLBCX, HLBCY, KRR, KSV,PDK2U, PDK4U, & PDK2TH, PDK4TH,PDK2SV, PDK4SV, KMI, & PUT, PVT, PWT, PTHT, PTKET, PRT, PSVT, & PLSUM,PLSVM,PLSWM,PLSTHM,PLSRVM,PRHODJ, & PRUS, PRVS, PRWS, PRTHS, PRTKES, PRRS, PRSVS, & OZDIFFU,ONUMDIFU, ONUMDIFTH, ONUMDIFSV, & TPHALO2LIST, TPHALO2LSLIST,PZDIFFU_HALO2, KTCOUNT ) ! ###################################################################### ! !!**** *NUM_DIFF2 * - routine to call the NUM_DIFF2_ALGO routine !! for each prognostic variable !! !! !! PURPOSE !! ------- !! The purpose of this routine is to call the NUM_DIFF2_ALGO !! routine for each prognostic variable. The NUM_DIFF2_ALGO routine !! applies an horizontal diffusion to the prognostic variables. !! !! !!** METHOD !! ------ !! For each prognostic variable the NUM_DIFF2 routine calls !! the NUM_DIFF2_ALGO routine which computes the numerical diffusion !! for a 3D field. !! The following variables are passed as argument to NUM_DIFF2_ALGO : !! !! -- The source term and the variable at t-dt !! !! -- The density PRHODJ !! or the mean operator in x or y or z direction applied to PRHODJ !! (as appropriate) !! !! -- The second layer of the halo of the field at t-dt !! !! -- The second layer of the halo of the LS field at t-dt, if necessary !! !! -- The LS field at t-dt, if necessary !! !! -- The localisation on the model grid : !! !! IGRID = 1 for mass grid point !! IGRID = 2 for U grid point !! IGRID = 3 for V grid point !! IGRID = 4 for w grid point !! !! EXTERNAL !! -------- !! BUDGET : Stores the different budget components !! (not used in current version) !! !! NUM_DIFF2_ALGO : applies an horizontal diffusion to the prognostic !! variables !! !! IMPLICIT ARGUMENTS !! ------------------ !! MODULE MODD_BUDGET: !! NBUMOD : model in which budget is calculated !! CBUTYPE : type of desired budget !! 'CART' for cartesian box configuration !! 'MASK' for budget zone defined by a mask !! 'NONE' ' for no budget !! NBUPROCCTR : process counter used for each budget variable !! Switches for budgets activations: !! !! LBU_RU : logical for budget of RU (wind component along x) !! !! LBU_RU : logical for budget of RU (wind component along x) !! .TRUE. = budget of RU !! .FALSE. = no budget of RU !! LBU_RV : logical for budget of RV (wind component along y) !! .TRUE. = budget of RV !! .FALSE. = no budget of RV !! LBU_RW : logical for budget of RW (wind component along z) !! .TRUE. = budget of RW !! .FALSE. = no budget of RW !! LBU_RTH : logical for budget of RTH (potential temperature) !! .TRUE. = budget of RTH !! .FALSE. = no budget of RTH !! LBU_RTKE : logical for budget of RTKE (turbulent kinetic energy) !! .TRUE. = budget of RTKE !! .FALSE. = no budget of RTKE !! LBU_RRV : logical for budget of RRV (water vapor) !! .TRUE. = budget of RRV !! .FALSE. = no budget of RRV !! LBU_RRC : logical for budget of RRC (cloud water) !! .TRUE. = budget of RRC !! .FALSE. = no budget of RRC !! LBU_RRR : logical for budget of RRR (rain water) !! .TRUE. = budget of RRR !! .FALSE. = no budget of RRR !! LBU_RRI : logical for budget of RRI (ice) !! .TRUE. = budget of RRI !! .FALSE. = no budget of RRI !! LBU_RRS : logical for budget of RRS (snow) !! .TRUE. = budget of RRS !! .FALSE. = no budget of RRS !! LBU_RRG : logical for budget of RRG (graupel) !! .TRUE. = budget of RRG !! .FALSE. = no budget of RRG !! LBU_RRH : logical for budget of RRH (hail) !! .TRUE. = budget of RRH !! .FALSE. = no budget of RRH !! LBU_RSV : logical for budget of RSVx (scalar variable) !! .TRUE. = budget of RSVx !! .FALSE. = no budget of RSVx !! !! MODULE MODD_PARALLEL2 !! !! HALO2_ll : type for a set of four arrays representing !! the second layer of the halo !! !! MODULE MODD_ARGSLIST !! !! HALO2LIST_ll : type for a list of "HALO2_lls" !! !! REFERENCE !! --------- !! !! Book2 of documentation ( routine NUM_DIFF2 ) !! !! AUTHOR !! ------ !! E. Richard * Laboratoire d'Aerologie* !! J.-P. Pinty * Laboratoire d'Aerologie* !! !! MODIFICATIONS !! ------------- !! Original 15/10/94 !! J. Stein 08/12/94 change the variable to be diffused rho*J*phi => phi !! J. Stein 20/03/95 remove R from the historical variables + switch 2D !! + switch for TKE unused !! 01/04/95 (Ph. Hereil J. Nicolau) add the budget computation !! 16/10/95 (J. Stein) change the budget calls !! 19/12/96 (J.-P. Pinty) update the budget calls !! 09/06/98 (P. Kloos) restructure to run in parallel !! 06/11/02 (V. Masson) update the budget calls !! 03/11/04 (G. Zängl) add calls for truly horizontal diffusion !! 05/06 (C.Lac) Remove EPS !! 05/07 (C.Lac) Separation between variables !! 07/09 (C.Lac) Correction on budget calls !! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 !! !------------------------------------------------------------------------------- ! !* 0. DECLARATIONS ! ------------ ! USE MODE_POS USE MODE_ll USE MODE_IO_ll USE MODD_ARGSLIST_ll, ONLY : LIST_ll ! USE MODD_PARAMETERS USE MODD_CONF USE MODD_BUDGET ! USE MODI_SHUMAN USE MODI_GDIV USE MODI_CONTRAV USE MODI_BUDGET USE MODI_GRADIENT_UV USE MODI_GRADIENT_VW USE MODI_GRADIENT_UW USE MODI_GRADIENT_M USE MODI_GRADIENT_U USE MODI_GRADIENT_V USE MODI_GRADIENT_W USE MODD_FIELD_n ! USE MODE_TYPE_ZDIFFU USE MODD_IBM_PARAM_n USE MODI_IBM_WRITE USE MODD_REF_n, ONLY: XRHODJ,XRHODREF USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZZ USE MODD_DYN_n, ONLY: XTSTEP USE MODD_CTURB USE MODD_CST USE MODD_METRICS_n, ONLY: XDXX,XDYY,XDZZ,XDZX,XDZY USE MODD_LBC_n ! IMPLICIT NONE ! !* 0.1 Declarations of dummy arguments : ! ! CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! X direction LBC type CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! Y direction LBC type ! INTEGER, INTENT(IN) :: KRR ! Number of moist variables INTEGER, INTENT(IN) :: KSV ! Number of Scalar Variables INTEGER, INTENT(IN) :: KMI ! Model index INTEGER, INTENT(IN) :: KTCOUNT ! REAL, INTENT(IN) :: PDK2U ! 2nd order dif. coef. /dx2 REAL, INTENT(IN) :: PDK4U ! 4th order dif. coef. /dx4 ! for momentum REAL, INTENT(IN) :: PDK2TH ! for theta and mixing ratios REAL, INTENT(IN) :: PDK4TH ! REAL, INTENT(IN) :: PDK2SV ! for theta and mixing ratios REAL, INTENT(IN) :: PDK4SV ! ! REAL, DIMENSION(:,:,:), INTENT(IN) :: PLSUM,PLSVM ! Larger REAL, DIMENSION(:,:,:), INTENT(IN) :: PLSWM,PLSTHM ! Scales REAL, DIMENSION(:,:,:), INTENT(IN) :: PLSRVM ! Variables ! REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry density of ! ! anelastic reference state * Jacobian ! REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT, PWT , PTKET REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT , PSVT ! REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS, PRWS ! Source terms REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTHS, PRTKES ! REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS , PRSVS ! ! ! Fields for z-diffusion LOGICAL, INTENT(IN) :: OZDIFFU LOGICAL, INTENT(IN) :: ONUMDIFU LOGICAL, INTENT(IN) :: ONUMDIFTH LOGICAL, INTENT(IN) :: ONUMDIFSV TYPE(HALO2LIST_ll), POINTER :: TPHALO2LIST, TPHALO2LSLIST ! list for diffusion TYPE(TYPE_ZDIFFU_HALO2) :: PZDIFFU_HALO2 ! ! !* 0.2 Declarations of local variables : ! INTEGER :: JRR ! Loop index for moist variables INTEGER :: JSV ! Loop index for Scalar Variables INTEGER:: IIB,IJB,IKB ! Begining useful area in x,y directions INTEGER:: IIE,IJE,IKE ! End useful area in x,y directions INTEGER :: IIU,IJU,IKU ! LOGICAL :: GTKEALLOC ! true if TKE arrays are not zero-sized ! TYPE(HALO2LIST_ll), POINTER :: TZHALO2LIST, TZHALO2LSLIST ! INTEGER :: IGRID ! localisation on the model grid REAL,DIMENSION(:,:,:), ALLOCATABLE:: ZTMQ,ZTMR,ZTMS,ZTME REAL,DIMENSION(:,:,:), ALLOCATABLE:: ZA2U,ZA2V,ZA2W,ZA4U,ZA4V,ZA4W REAL,DIMENSION(:,:,:), ALLOCATABLE:: ZBU,ZBV,ZBW,ZTEU,ZTEV,ZTEW REAL :: ZCOEFA,ZCOEFB,ZCOEFC,ZCOEFD,ZIBM_CFL TYPE(LIST_ll), POINTER :: TZFIELDC_ll INTEGER :: IINFO_ll INTEGER :: IIE2,IJE2,IKE2 INTEGER :: IIE1,IJE1,IKE1,IIB1,IJB1,IKB1 ! !------------------------------------------------------------------------------- ! constants ZCOEFA = PDK4U ! ponderation of the artificial diffusion ZCOEFB = PDK2U ! ponderation of the presence of the turbulent diffusion ZCOEFC = 1. ! decrease of the diffusion near fluid-solid interface ZCOEFD = 2. ! increase of the diffusion in regard of the non-linearity ! !------------------------------------------------------------------------------- ! initialization CALL GET_INDICE_ll(IIB,IJB,IIE,IJE) IIU=SIZE(PUT,1) IJU=SIZE(PUT,2) IKU=SIZE(PUT,3) IKB=2 IKE=IKU-1 ! !------------------------------------------------------------------------------- ! beta_i computation ALLOCATE(ZBU(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) ALLOCATE(ZBV(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) ALLOCATE(ZBW(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) IIB1 = IIB-1 IJB1 = IJB-1 IF (LWEST_ll()) IIB1=IIB+1 IF (LSOUTH_ll()) IJB1=IJB+1 IKB1=IKB!+1 IIE1 = IIE+1 IJE1 = IJE+1 IF (LEAST_ll()) IIE1=IIE-1 IF (LNORTH_ll()) IJE1=IJE-1 IKE1=IKE-1 ZBU = 1. ZBV = 1. ZBW = 1. WHERE (ABS(PUT(IIB1+1:IIE1+1,:,:)-PUT(IIB1:IIE1,:,:)).GT.XIBM_EPSI) ZBU(IIB1:IIE1,:,:)=(1./3.)*(PUT(IIB1+2:IIE1+2,:,:)-PUT(IIB1-1:IIE1-1,:,:))/(PUT(IIB1+1:IIE1+1,:,:)-PUT(IIB1:IIE1,:,:)) WHERE (ABS(PVT(:,IJB1+1:IJE1+1,:)-PVT(:,IJB1:IJE1,:)).GT.XIBM_EPSI) ZBV(:,IJB1:IJE1,:)=(1./3.)*(PVT(:,IJB1+2:IJE1+2,:)-PVT(:,IJB1-1:IJE1-1,:))/(PVT(:,IJB1+1:IJE1+1,:)-PVT(:,IJB1:IJE1,:)) WHERE (ABS(PWT(:,:,IKB1+1:IKE1+1)-PWT(:,:,IKB1:IKE1)).GT.XIBM_EPSI) ZBW(:,:,IKB1:IKE1)=(1./3.)*(PWT(:,:,IKB1+2:IKE1+2)-PWT(:,:,IKB1-1:IKE1-1))/(PWT(:,:,IKB1+1:IKE1+1)-PWT(:,:,IKB1:IKE1)) ! !------------------------------------------------------------------------------- ! alpha_i computation (parameter to estimate the non-linearity) ZBU = MIN(ZBU,+4.-XIBM_EPSI) ZBV = MIN(ZBV,+4.-XIBM_EPSI) ZBW = MIN(ZBW,+4.-XIBM_EPSI) ZBU = (5./2.)*(1.-3./(4.-ZBU)) ZBV = (5./2.)*(1.-3./(4.-ZBV)) ZBW = (5./2.)*(1.-3./(4.-ZBW)) ZBU = MIN(ZBU,1.) ZBV = MIN(ZBV,1.) ZBW = MIN(ZBW,1.) ZBU = MAX(ZBU,0.) ZBV = MAX(ZBV,0.) ZBW = MAX(ZBW,0.) ! ZBW(:,:,IKB-1) = 0.!ZBW(:,:,IKB) ! ZBW(:,:,IKB) = 0.!ZBW(:,:,IKB-1) = 0.! ! ZBW(:,:,IKE+1) = 0.!ZBW(:,:,IKE) !IF (LSOUTH_ll()) ZBV(:,IJB-1,:) = 0.!ZBV(:,IJB,:) !IF (LSOUTH_ll()) ZBV(:,IJB ,:) = 0.!ZBV(:,IJB,:) !IF (LNORTH_ll()) ZBV(:,IJE+1,:) = 0.!ZBV(:,IJE,:) !IF ( LWEST_ll()) ZBU(IIB-1,:,:) = 0.!ZBU(IIB,:,:) !IF ( LWEST_ll()) ZBU(IIB ,:,:) = 0.!ZBU(IIB,:,:) !IF ( LEAST_ll()) ZBU(IIE+1,:,:) = 0.!ZBU(IIE,:,:) ! !------------------------------------------------------------------------------- ! location of the artificial term ALLOCATE(ZTMR(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) !ALLOCATE(ZTMQ(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) ZTMR(:,:,:)=1. !ZTMQ(:,:,:)=0. IIB1 = IIB-2 IJB1 = IJB-2 IF (LWEST_ll()) IIB1=IIB+1 IF (LSOUTH_ll()) IJB1=IJB+1 IKB1=IKB!+1 IIE1 = IIE+2 IJE1 = IJE+2 IF (LEAST_ll()) IIE1=IIE-1 IF (LNORTH_ll()) IJE1=IJE-1 IKE1=IKE-1 !ZTMR(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)= -(XIBM_LS(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1,1)/& ! (XRHODJ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1) /& ! XRHODREF(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1))**(1./3.)-1.)*ZCOEFC !ZTMR(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)= (0.5*(XZZ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)+ XZZ(IIB1:IIE1,IJB1:IJE1,IKB1+1:IKE1+1))/& ! (XRHODJ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)/& ! XRHODREF(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1))**(1./3.)-1.)*ZCOEFC !ZTMQ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)= 1.-exp(-( XIBM_LS(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1,1)/& ! (XRHODJ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1) /& ! XRHODREF(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1))**(1./3.)*ZCOEFC)**2.) !ZTMR(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)= 1.-exp(-(0.5*(XZZ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)+ XZZ(IIB1:IIE1,IJB1:IJE1,IKB1+1:IKE1+1))/& ! (XRHODJ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)/& ! XRHODREF(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1))**(1./3.)*ZCOEFC)**2.) !ZTMQ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)= MIN(ZTMQ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1),ZTMR(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)) !ZTMQ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)= MIN(ZTMQ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1), 1.) !ZTMQ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)= MAX(ZTMQ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1), 0.) !ZTMR(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)=1./2.*(1.-cos(3.14159*ZTMQ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1))) !ZTMR(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)= MIN(ZTMQ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1),ZTMR(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)) ! ZTMR(:,:,IKB-1) = 0. ! ZTMR(:,:,IKE ) = 1. ! ZTMR(:,:,IKE+1) = 1. !IF (LSOUTH_ll()) ZTMR(:,IJB-1,:) = 1. !IF (LNORTH_ll()) ZTMR(:,IJE+1,:) = 1. !IF ( LWEST_ll()) ZTMR(IIB-1,:,:) = 1. !IF ( LEAST_ll()) ZTMR(IIE+1,:,:) = 1. !DEALLOCATE(ZTMQ) ! !------------------------------------------------------------------------------- ! artificial viscosity/hyperviscosity minus molecular/turbulent viscosity (depending on the non-linearity) ALLOCATE(ZTME(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) ZTME(:,:,:) = ZCOEFB*(2.*PTKET(:,:,:))**0.5/(XRHODJ(:,:,:)/XRHODREF(:,:,:))**(1./3.)+& XIBM_VISC_U/(XRHODJ(:,:,:)/XRHODREF(:,:,:))**(2./3.) ! ZTME(:,:,IKB-1) = ZTME(:,:,IKB) ! ZTME(:,:,IKE+1) = ZTME(:,:,IKE) !IF (LSOUTH_ll()) ZTME(:,IJB-1,:) = ZTME(:,IJB,:) !IF (LNORTH_ll()) ZTME(:,IJE+1,:) = ZTME(:,IJE,:) !IF ( LWEST_ll()) ZTME(IIB-1,:,:) = ZTME(IIB,:,:) !IF ( LEAST_ll()) ZTME(IIE+1,:,:) = ZTME(IIE,:,:) ALLOCATE(ZTEU(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) ALLOCATE(ZTEV(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) ALLOCATE(ZTEW(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) ZIBM_CFL=MIN(1.,XIBM_CFL) ZTEU = MAX(0.,ZBU**ZCOEFD*ZCOEFA*ZIBM_CFL**2./(4.*XTSTEP)- MXM(ZTME)) ZTEV = MAX(0.,ZBV**ZCOEFD*ZCOEFA*ZIBM_CFL**2./(4.*XTSTEP)- MYM(ZTME)) ZTEW = MAX(0.,ZBW**ZCOEFD*ZCOEFA*ZIBM_CFL**2./(4.*XTSTEP)- MZM(1,IKU,1,ZTME)) ZTEU = MIN(ZTEU,1./4./XTSTEP)*MXM( ZTMR) ZTEV = MIN(ZTEV,1./4./XTSTEP)*MYM( ZTMR) ZTEW = MIN(ZTEW,1./4./XTSTEP)*MZM(1,IKU,1,ZTMR) DEALLOCATE(ZTME,ZTMR) ! !------------------------------------------------------------------------------- ! coeff if non regular cartesian mesh (z direction) ALLOCATE(ZTMS(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) ZTMS = 1. ZTMS(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1) = (XRHODJ(IIB1:IIE1,IJB1:IJE1,IKB1 :IKE1 ) /XRHODREF(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1))**(1./3.)/& ( XZZ(IIB1:IIE1,IJB1:IJE1,IKB1+1:IKE1+1) -XZZ(IIB1:IIE1,IJB1:IJE1,IKB1:IKE1)) ! ZTMS(:,:,IKB-1) = ZTMS(:,:,IKB) ! ZTMS(:,:,IKE+1) = ZTMS(:,:,IKE) !IF (LSOUTH_ll()) ZTMS(:,IJB-1,:) = ZTMS(:,IJB,:) !IF (LNORTH_ll()) ZTMS(:,IJE+1,:) = ZTMS(:,IJE,:) !IF ( LWEST_ll()) ZTMS(IIB-1,:,:) = ZTMS(IIB,:,:) !IF ( LEAST_ll()) ZTMS(IIE+1,:,:) = ZTMS(IIE,:,:) ! !------------------------------------------------------------------------------- ! coefficient for viscous and hyperviscous term ALLOCATE(ZA2U(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) ALLOCATE(ZA2V(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) ALLOCATE(ZA2W(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) ALLOCATE(ZA4U(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) ALLOCATE(ZA4V(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) ALLOCATE(ZA4W(SIZE(PRUS,1),SIZE(PRVS,2),SIZE(PRWS,3))) ZA4U = (1./4.)*(1.-ZBU)*ZTEU ZA4V = (1./4.)*(1.-ZBV)*ZTEV ZA4W = (1./4.)*(1.-ZBW)*ZTEW*ZTMS**4. ZA2U = (1./1.)*(0.+ZBU)*ZTEU+XIBM_VISC_U/(XRHODJ/XRHODREF)**(2./3.) ZA2V = (1./1.)*(0.+ZBV)*ZTEV+XIBM_VISC_U/(XRHODJ/XRHODREF)**(2./3.) ZA2W =((1./1.)*(0.+ZBW)*ZTEW+XIBM_VISC_U/(XRHODJ/XRHODREF)**(2./3.))*ZTMS**2. DEALLOCATE(ZTEU,ZTEV,ZTEW,ZBU,ZBV,ZBW) !IGRID = IKB !ZA4W(:,:,IGRID) = 0. !ZA2W(:,:,IGRID) = 0. !IF (LSOUTH_ll()) THEN !IGRID = IJB !ZA4V(:,IGRID,:) = 0. !ZA2V(:,IGRID,:) = 0. !ENDIF !IF (LNORTH_ll()) THEN !IGRID = IJE+1 !ZA4V(:,IGRID,:) = 0. !ZA2V(:,IGRID,:) = 0. !ENDIF !IF (LWEST_ll()) THEN !IGRID = IIB !ZA4U(IGRID,:,:) = 0. !ZA2U(IGRID,:,:) = 0. !ENDIF !IF (LEAST_ll()) THEN !IGRID = IIE+1 !ZA4U(IGRID,:,:) = 0. !ZA2U(IGRID,:,:) = 0. !ENDIF ! !------------------------------------------------------------------------------- ! compute the diffusion term IGRID = 2 CALL NUM_DIFF2_ALGO( PRUS, PUT, IGRID, XRHODJ2(:,:,:,1),MXM(ZA2U),MXM(ZA4U),MXM(ZA2V),MXM(ZA4V),MXM(ZA2W),MXM(ZA4W)) IGRID = 3 CALL NUM_DIFF2_ALGO( PRVS, PVT, IGRID, XRHODJ2(:,:,:,2),MYM(ZA2U),MYM(ZA4U),MYM(ZA2V),MYM(ZA4V),MYM(ZA2W),MYM(ZA4W)) IGRID = 4 CALL NUM_DIFF2_ALGO( PRWS, PWT, IGRID, XRHODJ2(:,:,:,3),MZM(1,IKU,1,ZA2U) ,MZM(1,IKU,1,ZA4U), MZM(1,IKU,1,ZA2V) ,& MZM(1,IKU,1,ZA4V) ,MZM(1,IKU,1,ZA2W), MZM(1,IKU,1,ZA4W) ) IF (LIBM.AND.LFLAT) THEN PRWS(:,:,IKB )=0. PRWS(:,:,IKB-1)=-PRWS(:,:,IKB+1) ENDIF ZA4U = 0. ZA4V = 0. ZA4W = 0. ZA2U(:,:,:) = XIBM_VISC_U/(XRHODJ(:,:,:)/XRHODREF(:,:,:))**(2./3.) ZA2V(:,:,:) = XIBM_VISC_U/(XRHODJ(:,:,:)/XRHODREF(:,:,:))**(2./3.) ZA2W(:,:,:) =(XIBM_VISC_U/(XRHODJ(:,:,:)/XRHODREF(:,:,:))**(2./3.))*ZTMS**2. !IGRID = IKB-1 !ZA2U(:,:,IGRID) = 0. !ZA2V(:,:,IGRID) = 0. !ZA2W(:,:,IGRID) = 0. !IGRID = IKE+1 !ZA2U(:,:,IGRID) = 0. !ZA2V(:,:,IGRID) = 0. !ZA2W(:,:,IGRID) = 0. !IF (LSOUTH_ll()) THEN !IGRID = IJB-1 !ZA2U(:,IGRID,:) = 0. !ZA2V(:,IGRID,:) = 0. !ZA2W(:,IGRID,:) = 0. !ENDIF !IF (LNORTH_ll()) THEN !IGRID = IJE+1 !ZA2U(:,IGRID,:) = 0. !ZA2V(:,IGRID,:) = 0. !ZA2W(:,IGRID,:) = 0 !ENDIF !IF (LWEST_ll()) THEN !IGRID = IIB-1 !ZA2U(IGRID,:,:) = 0. !ZA2V(IGRID,:,:) = 0. !ZA2W(IGRID,:,:) = 0. !ENDIF !IF (LEAST_ll()) THEN !IGRID = IIE+1 !ZA2U(IGRID,:,:) = 0. !ZA2V(IGRID,:,:) = 0. !ZA2W(IGRID,:,:) = 0. !ENDIF IGRID = 1 CALL NUM_DIFF2_ALGO(PRTKES, PTKET, IGRID, (PRHODJ), ZA2U, ZA4U, ZA2V, ZA4V, ZA2W, ZA4W) ZA2U(:,:,:) = XIBM_VISC_T/XCPD/(XRHODJ(:,:,:)/XRHODREF(:,:,:))**(2./3.) ZA2V(:,:,:) = XIBM_VISC_T/XCPD/(XRHODJ(:,:,:)/XRHODREF(:,:,:))**(2./3.) ZA2W(:,:,:) =(XIBM_VISC_T/XCPD/(XRHODJ(:,:,:)/XRHODREF(:,:,:))**(2./3.))*ZTMS**2. CALL NUM_DIFF2_ALGO(PRTHS , PTHT, IGRID, (PRHODJ), ZA2U, ZA4U, ZA2V, ZA4V, ZA2W, ZA4W) IF (LIBM_SOLAR) THEN XIBM_LAMT(:,:,1 ,1) = XIBM_LAMT(:,:, 2,1) XIBM_LAMT(:,:,IKU,1) = XIBM_LAMT(:,:,IKU-1,1) ZA2U(:,:,:) = 2. * MXM( XIBM_LAMT(:,:,:,1)) * GX_M_U(1,IKU,1,PTHT(:,:,:),XDXX,XDZZ,XDZX) ZA2V(:,:,:) = 2. * MYM( XIBM_LAMT(:,:,:,1)) * GY_M_V(1,IKU,1,PTHT(:,:,:),XDYY,XDZZ,XDZY) ZA2W(:,:,:) = 2. * MZM(1,IKU,1,XIBM_LAMT(:,:,:,1)) * GZ_M_W(1,IKU,1,PTHT(:,:,:),XDZZ) ZA2W(:,:,IKB) = 0. ZA2W(:,:,IKB-1) = 0. ZA4V(:,:,:) = 0. CALL GDIV(HLBCX,HLBCY,XDXX,XDYY,XDZX,XDZY,XDZZ,ZA2U(:,:,:),ZA2V(:,:,:),ZA2W(:,:,:),ZA4V(:,:,:)) WHERE (XIBM_LS(:,:,:,1).GT.-XIBM_EPSI) ZA4V(:,:,:) = 0. IF (LWEST_ll ()) ZA4V(1 ,:,:) = ZA4V(2 ,:,:) IF (LEAST_ll ()) ZA4V(IIU,:,:) = ZA4V(IIU-1,:,:) IF (LSOUTH_ll()) ZA4V(:,1 ,:) = ZA4V(:,2 ,:) IF (LNORTH_ll()) ZA4V(:,IJU,:) = ZA4V(:,IJU-1,:) ZA4V(:,:,1 ) = ZA4V(:,:, 2) ZA4V(:,:,IKU) = ZA4V(:,:,IKU-1) PRTHS(:,:,:) = PRTHS(:,:,:)+ZA4V(:,:,:)*(XRHODJ(:,:,:)/XRHODREF(:,:,:)) PRTHS(:,:,IKB-1) = PRTHS(:,:,IKB) XIBM_LAMT(:,:,:,1) = 0. ENDIF !IF (LIBM_SOLAR) THEN ! XIBM_LAMT(:,:,1 ,3) = XIBM_LAMT(:,:, 2,3) ! XIBM_LAMT(:,:,IKU,3) = XIBM_LAMT(:,:,IKU-1,3) ! ! ZA4U(:,:,:) = MXF(PUT) ! ZA2U(:,:,:) = 2. * MXM( XIBM_LAMT(:,:,:,3)) * GX_M_U(1,IKU,1,ZA4U(:,:,:),XDXX,XDZZ,XDZX) ! ZA2V(:,:,:) = 2. * MYM( XIBM_LAMT(:,:,:,3)) * GY_M_V(1,IKU,1,ZA4U(:,:,:),XDYY,XDZZ,XDZY) ! ZA2W(:,:,:) = 2. * MZM(1,IKU,1,XIBM_LAMT(:,:,:,3)) * GZ_M_W(1,IKU,1,ZA4U(:,:,:),XDZZ) ! ZA2W(:,:,IKB) = 0. ! ZA2W(:,:,IKB-1) = 0. ! ZA4V(:,:,:) = 0. ! CALL GDIV(HLBCX,HLBCY,XDXX,XDYY,XDZX,XDZY,XDZZ,ZA2U(:,:,:),ZA2V(:,:,:),ZA2W(:,:,:),ZA4V(:,:,:)) ! WHERE (XIBM_LS(:,:,:,1).GT.-XIBM_EPSI) ZA4V(:,:,:) = 0. ! IF (LWEST_ll ()) ZA4V(1 ,:,:) = ZA4V(2 ,:,:) ! IF (LEAST_ll ()) ZA4V(IIU,:,:) = ZA4V(IIU-1,:,:) ! IF (LSOUTH_ll()) ZA4V(:,1 ,:) = ZA4V(:,2 ,:) ! IF (LNORTH_ll()) ZA4V(:,IJU,:) = ZA4V(:,IJU-1,:) ! ZA4V(:,:,2 ) = 0. ! ZA4V(:,:,1 ) = ZA4V(:,:, 2) ! ZA4V(:,:,IKU) = ZA4V(:,:,IKU-1) ! PRUS(:,:,:) = PRUS(:,:,:)+MXM(ZA4V(:,:,:))*(XRHODJ(:,:,:)/XRHODREF(:,:,:)) ! PRUS(:,:,IKB-1) = PRUS(:,:,IKB) ! ! ZA4U(:,:,:) = MYF(PVT) ! ZA2U(:,:,:) = 2. * MXM( XIBM_LAMT(:,:,:,3)) * GX_M_U(1,IKU,1,ZA4U(:,:,:),XDXX,XDZZ,XDZX) ! ZA2V(:,:,:) = 2. * MYM( XIBM_LAMT(:,:,:,3)) * GY_M_V(1,IKU,1,ZA4U(:,:,:),XDYY,XDZZ,XDZY) ! ZA2W(:,:,:) = 2. * MZM(1,IKU,1,XIBM_LAMT(:,:,:,3)) * GZ_M_W(1,IKU,1,ZA4U(:,:,:),XDZZ) ! ZA2W(:,:,IKB) = 0. ! ZA2W(:,:,IKB-1) = 0. ! ZA4V(:,:,:) = 0. ! CALL GDIV(HLBCX,HLBCY,XDXX,XDYY,XDZX,XDZY,XDZZ,ZA2U(:,:,:),ZA2V(:,:,:),ZA2W(:,:,:),ZA4V(:,:,:)) ! WHERE (XIBM_LS(:,:,:,1).GT.-XIBM_EPSI) ZA4V(:,:,:) = 0. ! IF (LWEST_ll ()) ZA4V(1 ,:,:) = ZA4V(2 ,:,:) ! IF (LEAST_ll ()) ZA4V(IIU,:,:) = ZA4V(IIU-1,:,:) ! IF (LSOUTH_ll()) ZA4V(:,1 ,:) = ZA4V(:,2 ,:) ! IF (LNORTH_ll()) ZA4V(:,IJU,:) = ZA4V(:,IJU-1,:) ! ZA4V(:,:,2 ) = 0. ! ZA4V(:,:,1 ) = ZA4V(:,:, 2) ! ZA4V(:,:,IKU) = ZA4V(:,:,IKU-1) ! PRVS(:,:,:) = PRVS(:,:,:)+MYM(ZA4V(:,:,:))*(XRHODJ(:,:,:)/XRHODREF(:,:,:)) ! PRVS(:,:,IKB-1) = PRVS(:,:,IKB) ! ! ZA4U(:,:,:) = MZF(1,IKU,1,PWT) ! ZA2U(:,:,:) = 2. * MXM( XIBM_LAMT(:,:,:,3)) * GX_M_U(1,IKU,1,ZA4U(:,:,:),XDXX,XDZZ,XDZX) ! ZA2V(:,:,:) = 2. * MYM( XIBM_LAMT(:,:,:,3)) * GY_M_V(1,IKU,1,ZA4U(:,:,:),XDYY,XDZZ,XDZY) ! ZA2W(:,:,:) = 2. * MZM(1,IKU,1,XIBM_LAMT(:,:,:,3)) * GZ_M_W(1,IKU,1,ZA4U(:,:,:),XDZZ) ! ZA2W(:,:,IKB) = 0. ! ZA2W(:,:,IKB-1) = 0. ! ZA4V(:,:,:) = 0. ! CALL GDIV(HLBCX,HLBCY,XDXX,XDYY,XDZX,XDZY,XDZZ,ZA2U(:,:,:),ZA2V(:,:,:),ZA2W(:,:,:),ZA4V(:,:,:)) ! WHERE (XIBM_LS(:,:,:,1).GT.-XIBM_EPSI) ZA4V(:,:,:) = 0. ! IF (LWEST_ll ()) ZA4V(1 ,:,:) = ZA4V(2 ,:,:) ! IF (LEAST_ll ()) ZA4V(IIU,:,:) = ZA4V(IIU-1,:,:) ! IF (LSOUTH_ll()) ZA4V(:,1 ,:) = ZA4V(:,2 ,:) ! IF (LNORTH_ll()) ZA4V(:,IJU,:) = ZA4V(:,IJU-1,:) ! ZA4V(:,:,2 ) = 0. ! ZA4V(:,:,1 ) = ZA4V(:,:, 2) ! ZA4V(:,:,IKU) = ZA4V(:,:,IKU-1) ! PRWS(:,:,:) = PRWS(:,:,:)+MZM(1,IKU,1,ZA4V(:,:,:))*(XRHODJ(:,:,:)/XRHODREF(:,:,:)) ! PRWS(:,:,IKB ) = 0. ! PRWS(:,:,IKB-1) = PRWS(:,:,IKB) ! ! XIBM_LAMT(:,:,:,3) = 0. !ENDIF DEALLOCATE(ZA2U,ZA2V,ZA2W,ZA4U,ZA4V,ZA4W,ZTMS) ! !------------------------------------------------------------------------------- ! ! CONTAINS ! ! ! ! ############################################################ SUBROUTINE NUM_DIFF2_ALGO(PRFIELDS, PFIELDM, KGRID, PRHODJ,ZDK2U, ZDK4U, ZDK2V, ZDK4V, ZDK2W, ZDK4W) ! ############################################################ !! !!**** *NUM_DIFF2_ALGO * - routine to apply an horizontal diffusion to the !! prognostic variables !! !! PURPOSE !! ------- !! The purpose of this routine is to add a 2sd or 4th order horizontal !! diffusion term to an advected variable. !! !!** METHOD !! ------ !! In case of cyclic LBCs, the routine adds a numerical diffusion !! term obtained by applying an horizontal nabla4 operator to the !! prognostic variable on each grid level. In case of open LBCs, !! the nabla4 operator degenerates to a nabla2 operator on the first !! ring inside the computationnal domain. !! In the current version no budget computation is performed. !! The "halo2" (or the second layer of the halo) of the prognostic !! variable is passed as argument. The "halo2" of the LS field !! is also passed as argument, if necessary. !! !! IMPLICIT ARGUMENTS !! ------------------ !! !! Module MODD_PARALLEL2 !! !! HALO2_ll : type for a set of arrays representing !! the second layer of the halo !! !! REFERENCE !! --------- !! !! Book2 of documentation ( routine NUM_DIFF2 ) !! User Interface for the MesoNH Parallel Package !! !! AUTHOR !! ------ !! !! E. Richard * Laboratoire d'Aerologie* !! J.-P. Pinty * Laboratoire d'Aerologie* !! !! MODIFICATIONS !! ------------- !! Original 15/10/94 !! J. Stein 08/12/94 change the variable to be diffused rho*J*phi => phi !! J. Stein 20/03/95 remove R from the historical variables + switch 2D !! + switch for TKE unused !! 01/04/95 (Ph. Hereil J. Nicolau) add the budget computation !! 16/10/95 (J. Stein) change the budget calls !! 19/12/96 (J.-P. Pinty) update the budget calls !! 19/06/98 (P. Kloos) restructure to run in parallel !! 12/10 (J.Escobar) replace DX4 & DY4 by inline code for reproductibility !! 12/10 (J.Escobar) reorder all operator from - to + index for reproductibility !! !------------------------------------------------------------------------------- ! !* 0. DECLARATIONS ! ------------ ! USE MODE_ll ! USE MODD_CONF USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll USE MODI_NABLA4 ! IMPLICIT NONE ! !* 0.1 Declarations of dummy arguments : ! REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRFIELDS ! source term REAL, DIMENSION(:,:,:), INTENT(IN) :: PFIELDM ! variable at t-dt INTEGER, INTENT(IN) :: KGRID ! C grid localisation REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! RHODJ at KGRID REAL, DIMENSION(:,:,:), INTENT(IN) :: ZDK2U,ZDK2V,ZDK2W ! 2nd order dif. coef. /dx2 REAL, DIMENSION(:,:,:), INTENT(IN) :: ZDK4U,ZDK4V,ZDK4W ! 4th order dif. coef. /dx4 ! !* 0.2 Declarations of local variables : ! INTEGER:: IW,IW2,IE,IE2,IS,IS2,IN,IN2,IB,IT ! Coordinate of forth order diffusion area ! !------------------------------------------------------------------------------- ! !* 1. CALCULATE THE NUMERICAL DIFFUSION TERM IN THE X DIRECTION ! IW=IIB IE=IIE IW2=IW-1 IE2=IE+1 ! IF (KGRID==2) IE2=IE+1 IF (LWEST_ll()) IW2=IW+2 IF (LEAST_ll()) IE2=IE-2 ! PRFIELDS(IW2:IE2,:,:) = PRFIELDS(IW2:IE2,:,:) - PRHODJ(IW2:IE2,:,:) *( & ZDK4U(IW2:IE2,:,:)*( PFIELDM(IW2-2:IE2-2,:,:) + PFIELDM(IW2+2:IE2+2,:,:) & -4.*( PFIELDM(IW2-1:IE2-1,:,:) + PFIELDM(IW2+1:IE2+1,:,:) ) & +6.* PFIELDM(IW2:IE2,:,:) ) -& ZDK2U(IW2:IE2,:,:)*( PFIELDM(IW2-1:IE2-1,:,:)-2.*PFIELDM(IW2:IE2,:,:) & + PFIELDM(IW2+1:IE2+1,:,:) )) IF (LWEST_ll()) THEN ! IF (KGRID/=2) THEN PRFIELDS(IW,:,:) = PRFIELDS(IW,:,:) + PRHODJ(IW,:,:) *(& ZDK4U(IW,:,:)*( PFIELDM(IW,:,:) -2.*PFIELDM(IW,:,:) + PFIELDM(IW+1,:,:) ) +& ZDK2U(IW,:,:)*( PFIELDM(IW,:,:) -2.*PFIELDM(IW,:,:) + PFIELDM(IW+1,:,:) )) ENDIF PRFIELDS(IW+1,:,:) = PRFIELDS(IW+1,:,:) - PRHODJ(IW+1,:,:) *( & ZDK4U(IW+1,:,:)*( PFIELDM(IW,:,:) + PFIELDM(IW+3,:,:) & -4.*( PFIELDM(IW,:,:) + PFIELDM(IW+2,:,:) ) & +6.* PFIELDM(IW+1,:,:) )- & ZDK2U(IW+1,:,:)*( PFIELDM(IW,:,:) -2.*PFIELDM(IW+1,:,:) + & PFIELDM(IW+2,:,:) )) ENDIF ! IF (LEAST_ll()) THEN ! IF (KGRID/=2) THEN PRFIELDS(IE,:,:) = PRFIELDS(IE,:,:) + PRHODJ(IE,:,: ) *( & ZDK4U(IE,:,:)*( PFIELDM(IE-1,:,:)-2.*PFIELDM(IE,:,:) + PFIELDM(IE,:,:) ) +& ZDK2U(IE,:,:)*( PFIELDM(IE-1,:,:)-2.*PFIELDM(IE,:,:) + PFIELDM(IE,:,:) )) ELSE PRFIELDS(IE,:,:) = PRFIELDS(IE,:,:) - PRHODJ(IE,:,:) *( & ZDK4U(IE,:,:)*( PFIELDM(IE-2,:,:) + PFIELDM(IE ,:,:) & -4.*( PFIELDM(IE-1,:,:) + PFIELDM(IE ,:,:) ) & +6.* PFIELDM(IE ,:,:) ) -& ZDK2U(IE,:,:)*( PFIELDM(IE-1,:,:) -2.*PFIELDM(IE,:,:) +& PFIELDM(IE,:,:) )) ENDIF PRFIELDS(IE-1,:,:) = PRFIELDS(IE-1,:,:) - PRHODJ(IE-1,:,:) *( & ZDK4U(IE-1,:,:)*( PFIELDM(IE-3,:,:) + PFIELDM(IE,:,:) & -4.*( PFIELDM(IE-2,:,:) + PFIELDM(IE,:,:) ) & +6.* PFIELDM(IE-1,:,:) ) -& ZDK2U(IE-1,:,:)*( PFIELDM(IE-2,:,:) -2.*PFIELDM(IE-1,:,:) +& PFIELDM(IE,:,:) )) ENDIF ! !------------------------------------------------------------------------------- ! !* 2. COMPUTES THE 4TH ORDER DIFFUSION TERMS IN THE Y DIRECTION: ! --------------------------------------------------------- ! ! IS=IJB IN=IJE IS2=IS-1 IN2=IN+1 ! IF (KGRID==3) IN2=IN+1 IF (LSOUTH_ll()) IS2=IS+2 IF (LNORTH_ll()) IN2=IN-2 ! PRFIELDS(:,IS2:IN2,:) = PRFIELDS(:,IS2:IN2,:) - PRHODJ(:,IS2:IN2,:) * ( & ZDK4V(:,IS2:IN2,:)*( PFIELDM(:,IS2-2:IN2-2,:) + PFIELDM(:,IS2+2:IN2+2,:) & -4.*( PFIELDM(:,IS2-1:IN2-1,:) + PFIELDM(:,IS2+1:IN2+1,:) ) & +6.* PFIELDM(:,IS2:IN2,:) ) -& ZDK2V(:,IS2:IN2,:)*( PFIELDM(:,IS2-1:IN2-1,:) -2.*PFIELDM(:,IS2:IN2,:) + & PFIELDM(:,IS2+1:IN2+1,:) )) IF (LSOUTH_ll()) THEN ! IF (KGRID/=3) THEN PRFIELDS(:,IS,:) = PRFIELDS(:,IS,:) + PRHODJ(:,IS,:) *( & ZDK4V(:,IS,:)*( PFIELDM(:,IS,:) -2.*PFIELDM(:,IS,:) + PFIELDM(:,IS+1,:) ) +& ZDK2V(:,IS,:)*( PFIELDM(:,IS,:) -2.*PFIELDM(:,IS,:) + PFIELDM(:,IS+1,:) )) ENDIF ! PRFIELDS(:,IS+1,:) = PRFIELDS(:,IS+1,:) - PRHODJ(:,IS+1,:) *( & ZDK4V(:,IS+1,:)*( PFIELDM(:,IS,:) + PFIELDM(:,IS+3,:) & -4.*( PFIELDM(:,IS,:) + PFIELDM(:,IS+2,:) ) & +6.* PFIELDM(:,IS+1,:) ) -& ZDK2V(:,IS+1,:)*( PFIELDM(:,IS,:) -2.*PFIELDM(:,IS+1,:) +& PFIELDM(:,IS+2,:) )) ENDIF ! IF (LNORTH_ll()) THEN ! IF (KGRID/=3) THEN PRFIELDS(:,IN,:) = PRFIELDS(:,IN,:) + PRHODJ(:,IN,: ) *( & ZDK4V(:,IN,:)*( PFIELDM(:,IN-1,:) -2.*PFIELDM(:,IN,:) + PFIELDM(:,IN,:) ) +& ZDK2V(:,IN,:)*( PFIELDM(:,IN-1,:) -2.*PFIELDM(:,IN,:) + PFIELDM(:,IN,:) )) ELSE PRFIELDS(:,IN ,:) = PRFIELDS(:,IN ,:) - PRHODJ(:,IN ,:) * ( & ZDK4V(:,IN,:)*( PFIELDM(:,IN-2,:) + PFIELDM(:,IN ,:) & -4.*( PFIELDM(:,IN-1,:) + PFIELDM(:,IN ,:) ) & +6.* PFIELDM(:,IN ,:) ) -& ZDK2V(:,IN,:)*( PFIELDM(:,IN-1,:) -2.*PFIELDM(:,IN,:) + & PFIELDM(:,IN,:) )) ENDIF ! PRFIELDS(:,IN-1,:) = PRFIELDS(:,IN-1,:) - PRHODJ(:,IN-1,:) * ( & ZDK4V(:,IN-1,:)*( PFIELDM(:,IN-3,:) + PFIELDM(:,IN,:) & -4.*( PFIELDM(:,IN-2,:) + PFIELDM(:,IN,:) ) & +6.* PFIELDM(:,IN-1,:) ) -& ZDK2V(:,IN-1,:)*( PFIELDM(:,IN-2,:) -2.*PFIELDM(:,IN-1,:) +& PFIELDM(:,IN,:) )) ENDIF ! !* 2. COMPUTES THE 4TH ORDER DIFFUSION TERMS IN THE Z DIRECTION: ! --------------------------------------------------------- ! IB=IKB IT=IKE ! PRFIELDS(:,:,IB+2:IT-2) = PRFIELDS(:,:,IB+2:IT-2) - PRHODJ(:,:,IB+2:IT-2) * ( & ZDK4W(:,:,IB+2:IT-2)*( PFIELDM(:,:,IB :IT-4)+ PFIELDM(:,:,IB+4:IT ) & -4.*( PFIELDM(:,:,IB+1:IT-3)+ PFIELDM(:,:,IB+3:IT-1) ) & +6.* PFIELDM(:,:,IB+2:IT-2) ) -& ZDK2W(:,:,IB+2:IT-2)*( PFIELDM(:,:,IB+1:IT-3)-2.*PFIELDM(:,:,IB+2:IT-2) +& PFIELDM(:,:,IB+3:IT-1) )) IF (KGRID/=4) THEN PRFIELDS(:,:,IB) = PRFIELDS(:,:,IB) + PRHODJ(:,:,IB) *( & ZDK4W(:,:,IB)*( PFIELDM(:,:,IB) -2.*PFIELDM(:,:,IB) + PFIELDM(:,:,IB+1) ) +& ZDK2W(:,:,IB)* ( PFIELDM(:,:,IB) -2.*PFIELDM(:,:,IB) + PFIELDM(:,:,IB+1) )) ENDIF ! PRFIELDS(:,:,IB+1) = PRFIELDS(:,:,IB+1) - PRHODJ(:,:,IB+1) * ( & ZDK4W(:,:,IB+1)*( PFIELDM(:,:,IB) + PFIELDM(:,:,IB+3) & -4.*( PFIELDM(:,:,IB) + PFIELDM(:,:,IB+2) ) & +6.* PFIELDM(:,:,IB+1) ) -& ZDK2W(:,:,IB+1)*( PFIELDM(:,:,IB+2)-2.*PFIELDM(:,:,IB+1) +& PFIELDM(:,:,IB) )) IF (KGRID/=4) THEN PRFIELDS(:,:,IT) = PRFIELDS(:,:,IT) + PRHODJ(:,:,IT) *( & ZDK4W(:,:,IT)*( PFIELDM(:,:,IT-1) -2.*PFIELDM(:,:,IT) + PFIELDM(:,:,IT) ) +& ZDK2W(:,:,IT)*( PFIELDM(:,:,IT-1) -2.*PFIELDM(:,:,IT) + PFIELDM(:,:,IT) )) ELSE PRFIELDS(:,:,IT) = PRFIELDS(:,:,IT ) - PRHODJ(:,:,IT ) * ( & ZDK4W(:,:,IT)*( PFIELDM(:,:,IT-2) + PFIELDM(:,:,IT) & -4.*( PFIELDM(:,:,IT-1) + PFIELDM(:,:,IT) ) & +6.* PFIELDM(:,:,IT ) ) -& ZDK2W(:,:,IT)*( PFIELDM(:,:,IT-1) -2.*PFIELDM(:,:,IT) + & PFIELDM(:,:,IT) )) ENDIF PRFIELDS(:,:,IT-1) = PRFIELDS(:,:,IT-1) - PRHODJ(:,:,IT-1) * ( & ZDK4W(:,:,IT-1)*( PFIELDM(:,:,IT-3) + PFIELDM(:,:,IT) & -4.*( PFIELDM(:,:,IT-2) + PFIELDM(:,:,IT) ) & +6.* PFIELDM(:,:,IT-1) ) -& ZDK2W(:,:,IT-1)*( PFIELDM(:,:,IT-2) -2.*PFIELDM(:,:,IT-1) + PFIELDM(:,:,IT) )) ! !------------------------------------------------------------------------------- END SUBROUTINE NUM_DIFF2_ALGO ! END SUBROUTINE NUM_DIFF2