!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 MODI_IBM_SOLAR_SHADOW ! ############################ ! INTERFACE ! SUBROUTINE IBM_SOLAR_SHADOW(P_ANGLEX,P_ANGLEY,P_ANGLEZ,PMAX,PEXP,PERR,PPHI,PRAD_LW,PRAD_SW,IBAND,PALB,PEMI,PTRAD) ! REAL ,INTENT(IN) :: P_ANGLEX,P_ANGLEY,P_ANGLEZ REAL ,INTENT(IN) :: PMAX,PEXP,PERR REAL, DIMENSION(:,:,:) ,INTENT(IN) :: PPHI REAL, DIMENSION(:,:) ,INTENT(INOUT) :: PRAD_LW REAL, DIMENSION(:,:,:) ,INTENT(INOUT) :: PRAD_SW INTEGER ,INTENT(IN) :: IBAND REAL, DIMENSION(:,:,:) ,INTENT(IN) :: PALB REAL, DIMENSION(:,:) ,INTENT(IN) :: PEMI REAL, DIMENSION(:,:) ,INTENT(IN) :: PTRAD ! !------------------------------------------------------------------------------ ! END SUBROUTINE IBM_SOLAR_SHADOW ! END INTERFACE ! END MODULE MODI_IBM_SOLAR_SHADOW ! ! ########################################################################################################### SUBROUTINE IBM_SOLAR_SHADOW(P_ANGLEX,P_ANGLEY,P_ANGLEZ,PMAX,PEXP,PERR,PPHI,PRAD_LW,PRAD_SW,IBAND,PALB,PEMI,PTRAD) ! ########################################################################################################### ! !! !!**** IBM_SOLAR_SHADOW searchs the shadows region in presence of IBM !! !! PURPOSE !! ------- ! ! Redistribution of the SW_down radiation: ! 1 - from shadow ground energies to enlightened walls ! 2 - from feet of walls to enlightened walls ! !! METHOD !! ------ !! ! SOLAR_SHADOW = 1.0 if cell is an enlightened region ! SOLAR_SHADOW = 0.0 if cell is a solid region (IBM) ! SOLAR_SHADOW = 0.5 if cell is a shadow region ! SOLAR_SHADOW is computed from the top of the domain and function of the ! angle solar. If IB is detected and function of the angle solar, ! adjacent cells are considered as in a shadow region. The information is ! transmitted to lower levels. ! ! LW_down is applied to the energy balance of the roofs, taking into account sigma.T^4 ! LW_down = sigma.T^4 is considered on IB walls ! ! PRAD_SW is related to the ground/roof (if presence) surface SW radiative fluxes ! ZRAD_SW is related to the ground with GCT to know the information of surface SW radiative fluxes at ghosts ! !! EXTERNAL !! -------- !! SUBROUTINE ? !! !! IMPLICIT ARGUMENTS !! ------------------ !! MODD_? !! !! REFERENCE !! --------- !! !! AUTHOR !! ------ !! Franck Auguste (CERFACS-AE) !! !! MODIFICATIONS !! ------------- !! Original 01/01/2018 !! !------------------------------------------------------------------------------ ! !**** 0. DECLARATIONS ! --------------- ! module USE MODE_POS USE MODE_ll USE MODE_IO_ll ! ! declaration USE MODD_IBM_PARAM_n USE MODD_PARAMETERS, ONLY: JPVEXT,JPHEXT USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZZ USE MODD_METRICS_n, ONLY: XDXX,XDYY,XDZZ,XDZX,XDZY USE MODD_ARGSLIST_ll, ONLY : LIST_ll ! ! interface USE MODI_SHUMAN USE MODI_IBM_WRITE USE MODI_IBM_AFFECTP_2D USE MODI_IBM_INTERPOS USE MODI_GRADIENT_M USE MODI_GRADIENT_U USE MODI_GRADIENT_V USE MODI_GRADIENT_W USE MODD_DYN_n, ONLY: XTSTEP USE MODD_CST ! IMPLICIT NONE ! !------------------------------------------------------------------------------ ! ! 0.1 declarations of arguments ! REAL ,INTENT(IN) :: P_ANGLEX,P_ANGLEY,P_ANGLEZ ! Solar angle in the referent frame REAL ,INTENT(IN) :: PMAX,PEXP,PERR ! Altitude max of potential shadow REAL, DIMENSION(:,:,:) ,INTENT(IN) :: PPHI ! Levelset function REAL, DIMENSION(:,:) ,INTENT(INOUT) :: PRAD_LW ! Tendency due to LW REAL, DIMENSION(:,:,:) ,INTENT(INOUT) :: PRAD_SW ! Tendency due to SW INTEGER ,INTENT(IN) :: IBAND ! Spectral band number in SW REAL, DIMENSION(:,:,:) ,INTENT(IN) :: PALB ! Albedo per spectral band REAL, DIMENSION(:,:) ,INTENT(IN) :: PEMI ! Emissivity REAL, DIMENSION(:,:) ,INTENT(IN) :: PTRAD ! !------------------------------------------------------------------------------ ! ! 0.2 declaration of local variables ! INTEGER :: IIB,IJB,IKB,IIE,IJE,IKE ! Local domain size INTEGER :: IIU,IJU,IKU ! Global domain size INTEGER :: JI,JJ,JK,JL,JMMM,JMM,JM,JI2,JJ2,JK2,JS ! Loop index TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! List of fields to exchange INTEGER :: IINFO_ll INTEGER :: JIMIN,JIMAX,JJMIN,JJMAX,JKMAX ! Max. level where possible shadow INTEGER :: JICOEF,JJCOEF ! Loop direction REAL :: Z_NODEX,Z_NODEY,Z_NODEZ,Z_NODEX2,Z_NODEY2,Z_NODEZ2 ! Nodes location REAL :: ZCOEF,ZCOEFALL ! Diffusion coefficients REAL, DIMENSION(:,:,:) , ALLOCATABLE :: ZXREF,ZYREF,ZZREF REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZIBMVECN ! LS-grad compared to the solar angle REAL, DIMENSION(:,:,:) , ALLOCATABLE :: ZTMP_SHADOW,ZRAD_SW ! Temporary array 3D REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZIBM_SHADOW ! Temporary array 3D per spectral band REAL, DIMENSION(:) , ALLOCATABLE :: ZIBM_RAD_ALL ! Sum of the SW power at shadow surface per spectral band REAL, DIMENSION(:,:) , ALLOCATABLE :: ZIBM_SHA_SUR,ZIBM_SHA_WAL ! Weighting for redistribution of energy at wall feets REAL :: ZIBM_SHA_ALL INTEGER :: JKMAX1,JKMAX2 ! Vertical index of roofs ghosts ! !------------------------------------------------------------------------------ ! !**** 0. PRELIMINARIES ! ---------------- ! IIU=SIZE(PPHI,1) IJU=SIZE(PPHI,2) IKU=SIZE(PPHI,3) ! CALL GET_INDICE_ll(IIB,IJB,IIE,IJE) IKE = IKU - JPVEXT IKB = 1 + JPVEXT ! JKMAX = IKB DO JK=IKB+1,IKE-3 IF (XZZ(IIB,IJB,JK).LT.PMAX) JKMAX=JK+1 ENDDO ! JIMIN=IIB JIMAX=IIE JICOEF = +1 JJMIN=IJB JJMAX=IJE JJCOEF = +1 IF (P_ANGLEX.LT.XIBM_EPSI) THEN JIMIN=IIE JIMAX=IIB JICOEF = -1 ENDIF IF (P_ANGLEY.LT.XIBM_EPSI) THEN JJMIN=IJE JJMAX=IJB JJCOEF = -1 ENDIF ! ALLOCATE(ZXREF(SIZE(XIBM_LS,1),SIZE(XIBM_LS,2),SIZE(XIBM_LS,3))) ALLOCATE(ZYREF(SIZE(XIBM_LS,1),SIZE(XIBM_LS,2),SIZE(XIBM_LS,3))) ALLOCATE(ZZREF(SIZE(XIBM_LS,1),SIZE(XIBM_LS,2),SIZE(XIBM_LS,3))) CALL IBM_INTERPOS(ZXREF,ZYREF,ZZREF,'P') !------------------------------------------------------------------------------ ! !**** 1. DEFINITION OF SHADOW AND ENLIGHTENED REGIONS ! ----------------------------------------------- ! ALLOCATE(ZTMP_SHADOW(SIZE(XIBM_LS,1),SIZE(XIBM_LS,2),SIZE(XIBM_LS,3))) ZTMP_SHADOW (:,:, :)=1. WHERE (PPHI(:,:,:).GT.XIBM_EPSI) ZTMP_SHADOW (:,:, :)=0. ZTMP_SHADOW (:,:,IKB-1)=ZTMP_SHADOW (:,:, IKB) ! DO JK=JKMAX,IKB,-1 DO JJ=JJMIN,JJMAX,JJCOEF DO JI=JIMIN,JIMAX,JICOEF ZTMP_SHADOW(JI,JJ,JK) = 0. IF (PPHI(JI,JJ,JK).GT.XIBM_EPSI) GOTO 111 Z_NODEX = ZXREF(JI,JJ,JK)!(XXHAT(JI+1)+XXHAT(JI))/2. Z_NODEY = ZYREF(JI,JJ,JK)!(XYHAT(JJ+1)+XYHAT(JJ))/2. Z_NODEZ = ZZREF(JI,JJ,JK)!(XZZ(JI,JJ,JK+1)+XZZ(JI,JJ,JK))/2. ZCOEFALL = 0. DO JI2=JI-1,JI+1 DO JJ2=JJ-1,JJ+1 DO JK2=JK+0,JK+1 Z_NODEX2 = ZXREF(JI2,JJ2,JK2)!(XXHAT(JI2+1)+XXHAT(JI2))/2. Z_NODEY2 = ZYREF(JI2,JJ2,JK2)!(XYHAT(JJ2+1)+XYHAT(JJ2))/2. Z_NODEZ2 = ZZREF(JI2,JJ2,JK2)!(XZZ(JI2,JJ2,JK2+1)+XZZ(JI2,JJ2,JK2))/2. ZCOEF = MAX(0.,(P_ANGLEX*(Z_NODEX-Z_NODEX2)+P_ANGLEY*(Z_NODEY-Z_NODEY2)+P_ANGLEZ*(Z_NODEZ-Z_NODEZ2))**PEXP) ZCOEFALL = ZCOEFALL + ZCOEF ZTMP_SHADOW(JI,JJ,JK) = ZTMP_SHADOW(JI,JJ,JK)+ZCOEF*ZTMP_SHADOW(JI2,JJ2,JK2) ENDDO ENDDO ENDDO ZTMP_SHADOW(JI,JJ,JK) = ZTMP_SHADOW(JI,JJ,JK)/ZCOEFALL ZTMP_SHADOW(JI,JJ,JK) = max(0.,ZTMP_SHADOW(JI,JJ,JK)) ZTMP_SHADOW(JI,JJ,JK) = min(1.,ZTMP_SHADOW(JI,JJ,JK)) IF (ZTMP_SHADOW(JI,JJ,JK).GT.(1.-PERR)) THEN ZTMP_SHADOW(JI,JJ,JK)=1.0 ELSEIF (ZTMP_SHADOW(JI,JJ,JK).LT.PERR) THEN ZTMP_SHADOW(JI,JJ,JK)=0.0 ELSE ZTMP_SHADOW(JI,JJ,JK)=0.5 ENDIF IF (JK==IKB) ZTMP_SHADOW (JI,JJ,JK-1)=ZTMP_SHADOW (JI,JJ,JK) 111 CONTINUE ENDDO ENDDO NULLIFY(TZFIELDS_ll) CALL ADD2DFIELD_ll(TZFIELDS_ll,ZTMP_SHADOW(:,:,JK)) CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll) CALL CLEANLIST_ll(TZFIELDS_ll) ENDDO DEALLOCATE(ZXREF,ZYREF,ZZREF) ! !**** 2. SW-CONTRIBUTION IN THE SHADOW REGION ! --------------------------------------- ! ALLOCATE(ZIBM_SHADOW(SIZE(XIBM_LS,1),SIZE(XIBM_LS,2),SIZE(XIBM_LS,3),IBAND)) ALLOCATE(ZIBM_RAD_ALL(IBAND)) ZIBM_SHADOW(:,:,:,:)=0. ZIBM_RAD_ALL(:)=0. DO JS=1,IBAND WHERE (ABS(ZTMP_SHADOW(:,:,2)-0.5).LT.XIBM_EPSI**0.5) ZIBM_SHADOW(:,:,2,JS)=PRAD_SW(:,:,JS)*(1.-PALB(JI,JJ,JS))*& (1.-8.*PEMI(JI,JJ)*XTSTEP/XIBM_SOLAR_CPG/(XIBM_SOLAR_DZG+XIBM_SOLAR_DZR)*XSTEFAN*PTRAD(JI,JJ)**3.) ! Cp building = Cp earth PRAD_SW(:,:,JS)=0. END WHERE ZIBM_RAD_ALL(JS) = SUM3D_ll(ZIBM_SHADOW(:,:,:,JS),IINFO_ll)*NIBM_LAYER_P*1. ENDDO ! !------------------------------------------------------------------------------ ! !**** 3. SW-CONTRIBUTION TO ENLIGHTENED WALLS ! --------------------------------------- ! ! LS-vec(n) computation ALLOCATE(ZIBMVECN(SIZE(XIBM_LS,1),SIZE(XIBM_LS,2),SIZE(XIBM_LS,3),4)) ZIBMVECN(:,:,:,1) = GX_U_M(1,IKU,1,XIBM_LS(:,:,:,2),XDXX,XDZZ,XDZX) ZIBMVECN(:,:,:,2) = GY_V_M(1,IKU,1,XIBM_LS(:,:,:,3),XDYY,XDZZ,XDZY) ZIBMVECN(:,:,:,3) = GZ_W_M(1,IKU,1,XIBM_LS(:,:,:,4),XDZZ) ZIBMVECN(:,:,:,4) = ZIBMVECN(:,:,:,1)**2.+ZIBMVECN(:,:,:,2)**2.+ZIBMVECN(:,:,:,3)**2. ZIBMVECN(:,:,:,4) = ZIBMVECN(:,:,:,4)**0.5 WHERE (ZIBMVECN(:,:,:,4).GT.XIBM_EPSI) ZIBMVECN(:,:,:,1) = ZIBMVECN(:,:,:,1) / ZIBMVECN(:,:,:,4) ZIBMVECN(:,:,:,2) = ZIBMVECN(:,:,:,2) / ZIBMVECN(:,:,:,4) ZIBMVECN(:,:,:,3) = ZIBMVECN(:,:,:,3) / ZIBMVECN(:,:,:,4) ENDWHERE ! Estimate of the weighting function vec(n) and solar angle ALLOCATE(ZIBM_SHA_SUR(SIZE(XIBM_LS,1),SIZE(XIBM_LS,2))) ALLOCATE(ZIBM_SHA_WAL(SIZE(XIBM_LS,1),SIZE(XIBM_LS,2))) ZTMP_SHADOW(:,:,:)= 0. ZIBM_SHA_SUR(:,:) = MAX(0.,-XIBM_SOLAR_ANZ) ZIBM_SHA_WAL(:,:) = MAX(0.,-XIBM_SOLAR_ANZ) DO JMM=1,NIBM_LAYER_P ! JM = size(NIBM_GHOST_P,1) JI = 0 JJ = 0 JK = 0 DO WHILE ((JI==0.and.JJ==0.and.JK==0).and.JM>0) JI = NIBM_GHOST_P(JM,JMM,1,1) JJ = NIBM_GHOST_P(JM,JMM,1,2) JK = NIBM_GHOST_P(JM,JMM,1,3) IF (JI==0.and.JJ==0.and.JK==0) JM = JM - 1 ENDDO JMMM = JM ! IF (JMMM<=0) GO TO 666 DO JM = 1,JMMM JI = NIBM_GHOST_P(JM,JMM,1,1) JJ = NIBM_GHOST_P(JM,JMM,1,2) JK = NIBM_GHOST_P(JM,JMM,1,3) IF (JI==0.or.JJ==0.or.JK==0) GO TO 667 ZTMP_SHADOW(JI,JJ,JK)=(XIBM_SOLAR_ANX*ZIBMVECN(JI,JJ,JK,1)+& XIBM_SOLAR_ANY*ZIBMVECN(JI,JJ,JK,2)+& XIBM_SOLAR_ANZ*ZIBMVECN(JI,JJ,JK,3))*(1.+ZIBMVECN(JI,JJ,JK,3)**3.) ZTMP_SHADOW(JI,JJ,JK)=MAX(0.,ZTMP_SHADOW(JI,JJ,JK)) ZIBM_SHA_WAL(JI,JJ) = ZIBM_SHA_WAL(JI,JJ)+ ZTMP_SHADOW(JI,JJ,JK) 667 CONTINUE ENDDO ENDDO 666 CONTINUE ZIBM_SHA_ALL = SUM3D_ll( ZTMP_SHADOW(:,:,:),IINFO_ll) ALLOCATE(ZRAD_SW(SIZE(XIBM_LS,1),SIZE(XIBM_LS,2),IBAND)) DO JS=1,IBAND ZRAD_SW(:,:,JS)=PRAD_SW(:,:,JS)*(1.-8.*PEMI(:,:)*XTSTEP/XIBM_SOLAR_CPG/(XIBM_SOLAR_DZG+XIBM_SOLAR_DZR)*XSTEFAN*PTRAD(:,:)**3.)*& (1.-ZIBM_SHA_SUR(:,:))/ZIBM_SHA_WAL(:,:) CALL IBM_AFFECTP_2D(ZRAD_SW(:,:,JS),NIBM_LAYER_P,2.,1.,'CL2','LAI','NEU','CST','CN3',0.) ENDDO ! redistribution ZIBM_SHADOW(:,:,:,:)=0. XIBM_SW0(:,:,:)=0. XIBM_LW0(:,:,:)=0. DO JMM=1,NIBM_LAYER_P ! JM = size(NIBM_GHOST_P,1) JI = 0 JJ = 0 JK = 0 DO WHILE ((JI==0.and.JJ==0.and.JK==0).and.JM>0) JI = NIBM_GHOST_P(JM,JMM,1,1) JJ = NIBM_GHOST_P(JM,JMM,1,2) JK = NIBM_GHOST_P(JM,JMM,1,3) IF (JI==0.and.JJ==0.and.JK==0) JM = JM - 1 ENDDO JMMM = JM ! IF (JMMM<=0) GO TO 777 DO JM = 1,JMMM JI = NIBM_GHOST_P(JM,JMM,1,1) JJ = NIBM_GHOST_P(JM,JMM,1,2) JK = NIBM_GHOST_P(JM,JMM,1,3) IF (JI==0.or.JJ==0.or.JK==0) GO TO 778 DO JS=1,IBAND ZIBM_SHADOW(JI,JJ,JK,JS) = (ZIBM_RAD_ALL(JS)/ZIBM_SHA_ALL+ZRAD_SW(JI,JJ,JS)/ZIBM_SHA_WAL(JI,JJ))*ZTMP_SHADOW(JI,JJ,JK) XIBM_SW0(JI,JJ,JK) = (ZIBM_SHADOW(JI,JJ,JK,JS)+PRAD_SW(JI,JJ,JS)*MAX(0.,-ZIBMVECN(JI,JJ,JK,3)**3.))*(1.-PALB(JI,JJ,JS)) ENDDO XIBM_LW0(JI,JJ,JK)= (PRAD_LW(JI,JJ)-XSTEFAN*XIBM_SOLAR_P(JM,JMM,3)**4.)*MAX(0.,0.-ZIBMVECN(JI,JJ,JK,3)**3.)*PEMI(JI,JJ)& +XIBM_SOLAR_EMI*( 0.-XSTEFAN*XIBM_SOLAR_P(JM,JMM,3)**4.)*MAX(0.,1.+ZIBMVECN(JI,JJ,JK,3)**3.)*PEMI(JI,JJ) JI2 = NIBM_IMAGE_P(JM,JMM,1,1,1,1) JJ2 = NIBM_IMAGE_P(JM,JMM,1,1,1,2) PRAD_LW(JI2,JJ2)=PRAD_LW(JI2,JJ2)-XIBM_SOLAR_EMI*( 0.-XSTEFAN*XIBM_SOLAR_P(JM,JMM,3)**4.)*MAX(0.,1.+ZIBMVECN(JI,JJ,JK,3)**3.)*PEMI(JI,JJ) 778 CONTINUE ENDDO ENDDO 777 CONTINUE ! !------------------------------------------------------------------------------ ! !**** 3. SW-CHANGE IN SHADOW REGIONS AND AT WALL FEETS ! ------------------------------------------------ ! DO JS=1,IBAND PRAD_SW(:,:,JS)= ZIBM_SHA_SUR(:,:)/ZIBM_SHA_WAL(:,:)*PRAD_SW(:,:,JS) ENDDO ! !**** X. DEALLOCATIONS/CLOSES ! ----------------------- ! DEALLOCATE(ZTMP_SHADOW,ZIBM_SHADOW,ZIBM_SHA_SUR,ZIBM_SHA_WAL) DEALLOCATE(ZIBM_RAD_ALL,ZRAD_SW,ZIBMVECN) ! RETURN ! END SUBROUTINE IBM_SOLAR_SHADOW