main_rain_ice.F90

PROGRAM MAIN_RAIN_ICE

USE XRD_GETOPTIONS
USE GETDATA_RAIN_ICE_MOD
USE COMPUTE_DIFF
USE MODD_CONF
USE MODD_DIMPHYEX,   ONLY: DIMPHYEX_t
USE MODD_CST,        ONLY: CST_t, CST
USE MODD_RAIN_ICE_DESCR, ONLY : RAIN_ICE_DESCR
USE MODD_RAIN_ICE_PARAM, ONLY : RAIN_ICE_PARAM
USE MODD_PARAM_ICE,      ONLY: PARAM_ICE
USE MODI_RAIN_ICE
USE MODI_INI_CST
USE MODD_BUDGET!, ONLY: TBUCONF_ASSOCIATE, TBUDGETDATA, NBUDGET_RH, TBUCONF
USE STACK_MOD
USE OMP_LIB
USE YOMHOOK, ONLY : LHOOK, DR_HOOK
USE PARKIND1, ONLY : JPRB, JPIM


IMPLICIT NONE

INTEGER      :: KLON 
INTEGER      :: KLEV
INTEGER      :: KRR

REAL,    ALLOCATABLE, DIMENSION(:,:,:,:,:) :: PRS, PRS_OUT
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:,:) :: PFPR, PFPR_OUT
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:,:) :: PRT   
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PDZZ     
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PRHODJ  
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PRHODREF
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PEXNREF 
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PEXNREF2
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PPABSM  
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PHLC_HRC
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PHLC_HCF
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PHLI_HRI
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PHLI_HCF
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PTHT    
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PSIGS   
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PCLDFR  
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PTHS    
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PEVAP, PEVAP_OUT
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: PCIT, PCIT_OUT
REAL,    ALLOCATABLE, DIMENSION(:,:,:)     :: PSEA  
REAL,    ALLOCATABLE, DIMENSION(:,:,:)     :: PTOWN  
REAL,    ALLOCATABLE, DIMENSION(:,:,:)     :: PINPRR, PINPRR_OUT
REAL,    ALLOCATABLE, DIMENSION(:,:,:)     :: PINPRS, PINPRS_OUT
REAL,    ALLOCATABLE, DIMENSION(:,:,:)     :: PINPRG, PINPRG_OUT
REAL,    ALLOCATABLE, DIMENSION(:,:,:)     :: ZINDEP, ZINDEP_OUT
REAL,    ALLOCATABLE, DIMENSION(:,:,:,:)   :: ZRAINFR, ZRAINFR_OUT
REAL,    ALLOCATABLE, DIMENSION(:,:,:)     :: ZINPRC, ZINPRC_OUT
LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: LLMICRO 

INTEGER :: NPROMA, NGPBLKS, NFLEVG
INTEGER :: IBL, JLON, JLEV

TYPE(DIMPHYEX_t)         :: D, D0
CHARACTER (LEN=4)   :: CSUBG_AUCV_RC
CHARACTER (LEN=80)  :: CSUBG_AUCV_RI
LOGICAL             :: OSEDIC 
CHARACTER (LEN=4)   :: CSEDIM  
CHARACTER (LEN=4)   :: CMICRO  
REAL                :: PTSTEP 
LOGICAL :: OWARM 
LOGICAL :: OCND2 
LOGICAL :: LCRIAUTI
REAL :: ZCRIAUTI, ZT0CRIAUTI, ZCRIAUTC
TYPE(TBUDGETDATA), DIMENSION(NBUDGET_RH) :: YLBUDGET
LOGICAL                  :: LLCHECK
LOGICAL                  :: LLCHECKDIFF
LOGICAL                  :: LLDIFF
INTEGER                  :: IBLOCK1, IBLOCK2
INTEGER                  :: ISTSZ, JBLK1, JBLK2
INTEGER                  :: NTID, ITID
INTEGER                  :: JRR


INTEGER :: IPROMA
REAL, ALLOCATABLE :: PSTACK(:,:)
TYPE (STACK) :: YLSTACK

REAL(KIND=8) :: TS,TE
REAL(KIND=8) :: TSC, TEC, TSD, TED, ZTC, ZTD 
INTEGER :: ITIME, NTIME
INTEGER :: IRANK, ISIZE
LOGICAL :: LLVERBOSE, LLSTAT, LLBIND
REAL (KIND=JPRB) :: ZHOOK_HANDLE
CHARACTER(LEN=32) :: CLTEXT

CALL INITOPTIONS ()
NGPBLKS = 150
CALL GETOPTION ("--blocks", NGPBLKS)
NPROMA = 32
CALL GETOPTION ("--nproma", NPROMA)
NFLEVG = -1
CALL GETOPTION ("--nflevg", NFLEVG)
CALL GETOPTION ("--check",  LLCHECK)
CALL GETOPTION ("--checkdiff",  LLCHECKDIFF)
IBLOCK1 = 1
CALL GETOPTION ("--check-block-1", IBLOCK1)
IBLOCK2 = NGPBLKS
CALL GETOPTION ("--check-block-2", IBLOCK2)
CALL GETOPTION ("--stat", LLSTAT)
NTIME = 1
CALL GETOPTION ("--times", NTIME)
CALL GETOPTION ("--verbose", LLVERBOSE)
CALL GETOPTION ("--bind", LLBIND)
CALL CHECKOPTIONS ()

LLDIFF = .FALSE.

IRANK = 0
ISIZE = 1
IF (LLBIND) THEN
  CALL LINUX_BIND      (IRANK, ISIZE)
  CALL LINUX_BIND_DUMP (IRANK, ISIZE)
ENDIF

CALL GETDATA_RAIN_ICE (NPROMA, NGPBLKS, NFLEVG, LLMICRO, PEXNREF, PDZZ, PRHODJ, PRHODREF, &
&PEXNREF2, PPABSM, PCIT, PCLDFR, PHLC_HRC, PHLC_HCF, PHLI_HRI, PHLI_HCF, PTHT, PRT, PTHS, &
&PRS, PSIGS, PSEA, PTOWN, PCIT_OUT, PRS_OUT, ZINPRC, ZINPRC_OUT, PINPRR, PINPRR_OUT, PEVAP, PEVAP_OUT, &
&PINPRS, PINPRS_OUT, PINPRG, PINPRG_OUT, ZINDEP, ZINDEP_OUT, ZRAINFR, ZRAINFR_OUT, PFPR, PFPR_OUT, LLVERBOSE)

KLEV = SIZE (PRS, 3)
KRR  = SIZE (PRS, 4)

IF (LLVERBOSE) PRINT *, " KLEV = ", KLEV, " KRR = ", KRR

PRINT *, " NPROMA = ", NPROMA, " KLEV = ", KLEV, " NGPBLKS = ", NGPBLKS

CMICRO='ICE3'

PTSTEP = 25.0000000000000
KRR = 6
OSEDIC = .TRUE.
OCND2 = .FALSE.
CSEDIM = 'STAT'
CSUBG_AUCV_RC = 'PDF'
CSUBG_AUCV_RI = 'NONE'
OWARM = .TRUE.

LCRIAUTI=.TRUE.
ZCRIAUTI=0.2E-3
ZT0CRIAUTI=-5.
ZCRIAUTC=0.1E-2

CALL INIT_PHYEX (20, OWARM, CMICRO, CSEDIM, &
            & LCRIAUTI, ZCRIAUTI, ZT0CRIAUTI, ZCRIAUTC)
DO JRR=1, NBUDGET_RH
  YLBUDGET(JRR)%NBUDGET=JRR
ENDDO

D0%NIT  = NPROMA
D0%NIB  = 1
D0%NIE  = NPROMA
D0%NJT  = 1
D0%NJB  = 1
D0%NJE  = 1
D0%NIJT = D0%NIT * D0%NJT
D0%NIJB = 1
D0%NIJE = NPROMA
D0%NKL  = -1
D0%NKT  = KLEV
D0%NKA  = KLEV
D0%NKU  = 1
D0%NKB  = KLEV 
D0%NKE  = 1
D0%NKTB = 1
D0%NKTE = KLEV

ISTSZ = NPROMA * 20 * KLEV
ALLOCATE (PSTACK (ISTSZ, NGPBLKS))

TS = OMP_GET_WTIME ()

ZTD = 0.
ZTC = 0.

IF (LHOOK) CALL DR_HOOK ('MAIN',0,ZHOOK_HANDLE)

DO ITIME = 1, NTIME

  TSD = OMP_GET_WTIME ()

!directives pas a jour !$acc data &
!directives pas a jour !$acc      & copyin  (D0, CST, ICEP, NEB, KRR, HFRAC_ICE, HCONDENS, HLAMBDA3, HBUNAME, OSUBG_COND, OSIGMAS, OCND2, HSUBG_MF_PDF, PTSTEP, LMFCONV, &
!directives pas a jour !$acc      &          ZSIGQSAT, PRHODJ, PEXNREF, PRHODREF, PSIGS, PMFCONV, PPABSM, ZZZ, PCF_MF, PRC_MF, PRI_MF, ZRS, ZICE_CLD_WGT) &
!directives pas a jour !$acc      & copy    (PRS, PTHS), &
!directives pas a jour !$acc      & copyout (PSRCS, PCLDFR, PHLC_HRC, PHLC_HCF, PHLI_HRI, PHLI_HCF) &
!directives pas a jour !$acc      & create  (PSTACK) 

  TSC = OMP_GET_WTIME ()

#ifdef USE_OPENMP
!$OMP PARALLEL PRIVATE (D, ITID, JBLK1, JBLK2, IPROMA, ISIZE)
#endif

#ifdef _OPENACC
JBLK1 = 1 
JBLK2 = NGPBLKS
#endif

#ifdef USE_OPENMP
NTID = OMP_GET_MAX_THREADS ()
ITID = OMP_GET_THREAD_NUM ()
JBLK1 = 1 +  (NGPBLKS * (ITID+0)) / NTID
JBLK2 =      (NGPBLKS * (ITID+1)) / NTID


!PRINT *, ITID, JBLK1, JBLK2

#endif

!$acc parallel loop gang vector private (YLSTACK, IBL, JLON, D) collapse (2)

  DO IBL = JBLK1, JBLK2


#ifdef _OPENACC
  DO JLON = 1, NPROMA
    D = D0
    D%NIB = JLON
    D%NIE = JLON
    D%NIJB = JLON
    D%NIJE = JLON
#endif

#ifdef USE_OPENMP
    D = D0
#endif

#ifdef USE_STACK
    YLSTACK%L = LOC (PSTACK (1, IBL))
    YLSTACK%U = YLSTACK%L + ISTSZ * KIND (PSTACK)
#else
    YLSTACK%L = 0
    YLSTACK%U = 0
#endif

IPROMA=COUNT(LLMICRO(D%NIB:D%NIE,D%NJB:D%NJE,D%NKTB:D%NKTE,IBL))
ISIZE=IPROMA
CALL RAIN_ICE (D, CST, PARAM_ICE, RAIN_ICE_PARAM, &
             & RAIN_ICE_DESCR, TBUCONF, &
             & IPROMA, ISIZE, &
             & OSEDIC=OSEDIC, OCND2=OCND2, HSEDIM=CSEDIM, &
             & HSUBG_AUCV_RC=CSUBG_AUCV_RC, HSUBG_AUCV_RI=CSUBG_AUCV_RI,&
             & OWARM=OWARM, &
             & PTSTEP=2*PTSTEP, &
             & KRR=KRR, ODMICRO=LLMICRO(:,:,:,IBL), PEXN=PEXNREF(:,:,:,IBL),            &
             & PDZZ=PDZZ(:,:,:,IBL), PRHODJ=PRHODJ(:,:,:,IBL), PRHODREF=PRHODREF(:,:,:,IBL),PEXNREF=PEXNREF2(:,:,:,IBL),&
             & PPABST=PPABSM(:,:,:,IBL), PCIT=PCIT(:,:,:,IBL), PCLDFR=PCLDFR(:,:,:,IBL),  &
             & PHLC_HRC=PHLC_HRC(:,:,:,IBL), PHLC_HCF=PHLC_HCF(:,:,:,IBL), &
             & PHLI_HRI=PHLI_HRI(:,:,:,IBL), PHLI_HCF=PHLI_HCF(:,:,:,IBL), &
             & PTHT=PTHT,PRVT=PRT(:,:,:,1,IBL),PRCT=PRT(:,:,:,2,IBL), &
             & PRRT=PRT(:,:,:,3,IBL), &
             & PRIT=PRT(:,:,:,4,IBL), PRST=PRT(:,:,:,5,IBL), &
             & PRGT=PRT(:,:,:,6,IBL),       &
             & PTHS=PTHS(:,:,:,IBL), PRVS=PRS(:,:,:,1,IBL),PRCS=PRS(:,:,:,2,IBL),&
             & PRRS=PRS(:,:,:,3,IBL),&
             & PRIS=PRS(:,:,:,4,IBL),PRSS= PRS(:,:,:,5,IBL),PRGS= PRS(:,:,:,6,IBL),&
             & PINPRC=ZINPRC(:,:,IBL),PINPRR=PINPRR(:,:,IBL),PEVAP3D=PEVAP(:,:,:,IBL),&
             & PINPRS=PINPRS(:,:,IBL), PINPRG=PINPRG(:,:,IBL), PINDEP=ZINDEP(:,:,IBL), PRAINFR=ZRAINFR(:,:,:,IBL), &
             & PSIGS=PSIGS(:,:,:,IBL), &
             & TBUDGETS=YLBUDGET, KBUDGETS=SIZE(YLBUDGET), &
             & PSEA=PSEA, PTOWN=PTOWN, PFPR=PFPR(:,:,:,:,IBL))

#ifdef _OPENACC
    ENDDO
#endif

  ENDDO

#ifdef USE_OPENMP
!$OMP END PARALLEL
#endif

!$acc end parallel loop

  TEC = OMP_GET_WTIME ()

!$acc end data

  TED = OMP_GET_WTIME ()

  ZTC = ZTC + (TEC - TSC)
  ZTD = ZTD + (TED - TSD)

ENDDO

IF (LHOOK) CALL DR_HOOK ('MAIN',1,ZHOOK_HANDLE)

TE = OMP_GET_WTIME()

WRITE (*,'(A,F8.2,A)') 'elapsed time : ',TE-TS,' s'
WRITE (*,'(A,F8.4,A)') '          i.e. ',1000.*(TE-TS)/(NPROMA*NGPBLKS)/NTIME,' ms/gp'

PRINT *, " ZTD = ", ZTD, ZTD / REAL (NPROMA*NGPBLKS*NTIME)
PRINT *, " ZTC = ", ZTC, ZTC / REAL (NPROMA*NGPBLKS*NTIME)


IF (LLCHECK .OR. LLSTAT .OR. LLCHECKDIFF) THEN
  DO IBL = IBLOCK1, IBLOCK2
    PRINT *, " IBL = ", IBL
    DO JRR=1, KRR
      WRITE (CLTEXT, '("PRS JRR=",I3.3)') JRR
      CALL DIFF3 (CLTEXT,      PRS_OUT       (:,:,:,JRR,IBL), PRS      (:,:,:,JRR,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
      IF(JRR>=2) THEN
        WRITE (CLTEXT, '("PFPR JRR=",I3.3)') JRR
        CALL DIFF3 (CLTEXT,     PFPR_OUT      (:,:,:,JRR,IBL), PFPR     (:,:,:,JRR,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
      ENDIF
    ENDDO
    CALL DIFF3 ("PCIT",     PCIT_OUT      (:,:,:,IBL), PCIT     (:,:,:,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
    CALL DIFF2 ("ZINPRC",   ZINPRC_OUT    (:,:,IBL),   ZINPRC   (:,:,IBL)  , LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
    CALL DIFF2 ("PINPRRRS", PINPRR_OUT    (:,:,IBL),   PINPRR   (:,:,IBL)  , LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
    CALL DIFF3 ("PEVAP",    PEVAP_OUT     (:,:,:,IBL), PEVAP    (:,:,:,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
    CALL DIFF2 ("PINPRS",   PINPRS_OUT    (:,:,IBL),   PINPRS   (:,:,IBL)  , LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
    CALL DIFF2 ("PINPRG",   PINPRG_OUT    (:,:,IBL),   PINPRG   (:,:,IBL)  , LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
    CALL DIFF2 ("ZINDEP",   ZINDEP_OUT    (:,:,IBL),   ZINDEP   (:,:,IBL)  , LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
    CALL DIFF3 ("ZRAINFR",  ZRAINFR_OUT   (:,:,:,IBL), ZRAINFR  (:,:,:,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
  ENDDO
ENDIF

IF (LLCHECKDIFF) THEN
  IF (LLDIFF) THEN
    PRINT*, "THERE ARE DIFF SOMEWHERE"
  ELSE
    PRINT*, "THERE IS NO DIFF AT ALL"
  ENDIF
ENDIF

STOP

CONTAINS

SUBROUTINE INIT_PHYEX(KULOUT,LDWARM,CMICRO,CCSEDIM,LDCRIAUTI,&
                   PCRIAUTI,PT0CRIAUTI,PCRIAUTC)

USE MODD_RAIN_ICE_DESCR
USE MODD_RAIN_ICE_PARAM
USE MODD_PARAM_ICE
USE MODD_TURB_N, ONLY: TURB_GOTO_MODEL, CSUBG_MF_PDF

USE MODI_INI_RAIN_ICE

IMPLICIT NONE
! -----------------------------------------------------------------------
!     DUMMY INTEGER SCALARS
INTEGER, INTENT (IN) :: KULOUT
LOGICAL, INTENT (IN) :: LDWARM
CHARACTER(4), INTENT (IN) :: CMICRO 
CHARACTER(4), INTENT (IN) :: CCSEDIM
LOGICAL, INTENT (IN) :: LDCRIAUTI
REAL, INTENT (IN) :: PCRIAUTI
REAL, INTENT (IN) :: PT0CRIAUTI
REAL, INTENT (IN) :: PCRIAUTC
!-----------------------------------------------------------------------
!    LOCAL VARIABLES
REAL :: ZCRI0, ZTCRI0    
! -----------------------------------------------------------------------

CALL INI_CST
CALL TURB_GOTO_MODEL(1,1)
CALL PARAM_ICE_ASSOCIATE
CALL TBUCONF_ASSOCIATE
LBU_ENABLE=.FALSE.                                                                                                       
LBUDGET_U=.FALSE.
LBUDGET_V=.FALSE.
LBUDGET_W=.FALSE.
LBUDGET_TH=.FALSE.
LBUDGET_TKE=.FALSE.
LBUDGET_RV=.FALSE.
LBUDGET_RC=.FALSE.
LBUDGET_RR=.FALSE.
LBUDGET_RI=.FALSE.
LBUDGET_RS=.FALSE.
LBUDGET_RG=.FALSE.
LBUDGET_RH=.FALSE.
LBUDGET_SV=.FALSE.

!        1. Set implicit default values for MODD_PARAM_ICE
LWARM=LDWARM
CPRISTINE_ICE='PLAT'
CSEDIM=CCSEDIM
CSUBG_AUCV_RC='PDF'
CSUBG_AUCV_RI='NONE'
CSUBG_RC_RR_ACCR='NONE'
CSUBG_RR_EVAP='NONE'
CSUBG_PR_PDF='SIGM'
CSUBG_MF_PDF='TRIANGLE'
! Snow riming                                                                                                                       
CSNOWRIMING='M90 '
XFRACM90=0.1 ! Fraction used for the Murakami 1990 formulation
!                                                                                                                                   
LFEEDBACKT=.TRUE. ! When .TRUE. feed back on temperature is taken into account
LEVLIMIT=.TRUE.   ! When .TRUE. water vapour pressure is limited by saturation
LNULLWETG=.TRUE.  ! When .TRUE. graupel wet growth is activated with null rate (to allow water shedding)
LWETGPOST=.TRUE.  ! When .TRUE. graupel wet growth is activated with positive temperature (to allow water shedding)
LNULLWETH=.TRUE.  ! Same as LNULLWETG but for hail
LWETHPOST=.TRUE.  ! Same as LWETGPOST but for hail
LCONVHG=.TRUE. ! TRUE to allow the conversion from hail to graupel
LCRFLIMIT=.TRUE. !True to limit rain contact freezing to possible heat exchange
CFRAC_ICE_ADJUST='S' ! Ice/liquid partition rule to use in adjustment
CFRAC_ICE_SHALLOW_MF='S' ! Ice/liquid partition rule to use in shallow_mf
LSEDIM_AFTER=.FALSE. ! Sedimentation done after microphysics
XSPLIT_MAXCFL=0.8
LDEPOSC=.FALSE.  ! water deposition on vegetation
XVDEPOSC=0.02    ! deposition speed (2 cm.s-1)
!
!        2. Set implicit default values for MODD_RAIN_ICE_DESCR 
!                     et MODD_RAIN_ICE_PARAM
!
CALL INI_RAIN_ICE (KULOUT, CMICRO)
!update values from namparar
IF (LDCRIAUTI) THEN

  XCRIAUTI=PCRIAUTI
  XCRIAUTC=PCRIAUTC
  XT0CRIAUTI=PT0CRIAUTI
  !second point to determine 10**(aT+b) law
  ZTCRI0=-40.0
  ZCRI0=1.25E-6
  
  XBCRIAUTI=-( LOG10(XCRIAUTI) - LOG10(ZCRI0)*PT0CRIAUTI/ZTCRI0 )&
                   *ZTCRI0/(XT0CRIAUTI-ZTCRI0)
  XACRIAUTI=(LOG10(ZCRI0)-XBCRIAUTI)/ZTCRI0
  
ENDIF
! -----------------------------------------------------------------------

END SUBROUTINE INIT_PHYEX

END PROGRAM