source: trunk/LMDZ.TITAN/libf/muphytitan/mmp_gcm.f90 @ 3604

! Copyright (2017,2022-2023) Université de Reims Champagne-Ardenne
! Contributors : J. Burgalat (GSMA, URCA), B. de Batz de Trenquelléon (GSMA, URCA)
! email of the author : jeremie.burgalat@univ-reims.fr
!
! This software is a computer program whose purpose is to compute
! microphysics processes using a two-moments scheme.
!
! This library is governed by the CeCILL license under French law and
! abiding by the rules of distribution of free software.  You can  use,
! modify and/ or redistribute the software under the terms of the CeCILL
! license as circulated by CEA, CNRS and INRIA at the following URL
! "http://www.cecill.info".
!
! As a counterpart to the access to the source code and  rights to copy,
! modify and redistribute granted by the license, users are provided only
! with a limited warranty  and the software's author,  the holder of the
! economic rights,  and the successive licensors  have only  limited
! liability.
!
! In this respect, the user's attention is drawn to the risks associated
! with loading,  using,  modifying and/or developing or reproducing the
! software by the user in light of its specific status of free software,
! that may mean  that it is complicated to manipulate,  and  that  also
! therefore means  that it is reserved for developers  and  experienced
! professionals having in-depth computer knowledge. Users are therefore
! encouraged to load and test the software's suitability as regards their
! requirements in conditions enabling the security of their systems and/or
! data to be ensured and,  more generally, to use and operate it in the
! same conditions as regards security.
!
! The fact that you are presently reading this means that you have had
! knowledge of the CeCILL license and that you accept its terms.

!! file: mmp_gcm.f90
!! summary: YAMMS interfaces for the LMDZ GCM.
!! author: J. Burgalat
!! date: 2017,2022,2023
!! modifications: B. de Batz de Trenquelléon
MODULE MMP_GCM
  !! Interface to YAMMS for the LMDZ GCM.
  USE MMP_GLOBALS
  USE MM_LIB
  USE CFGPARSE
  USE DATASETS
  IMPLICIT NONE

  CONTAINS

  SUBROUTINE mmp_initialize(dt,p_prod,tx_prod,rc_prod,rplanet,g0, air_rad,air_mmol,clouds,cfgpath)
    !! Initialize global parameters of the model.
    !!
    !! The function initializes all the global parameters of the model from direct input.
    !! Boolean (and Fiadero) parameters are optional as they are rather testing parameters. Their
    !! default values are suitable for production runs.
    !! @note
    !! If the subroutine fails to initialize parameters, the run is aborted.
    !!
    !! @warning
    !! If OpenMP is activated, this subroutine must be called in an $OMP SINGLE statement as it
    !! initializes global variable that are not thread private.
    !!
    !! ```
    !! !$OMP SINGLE
    !! call mmp_initialize(...)
    !! !$OMP END SINGLE
    !! ```
    !!
    REAL(kind=mm_wp), INTENT(in)           :: dt
      !! Microphysics timestep in seconds.
    REAL(kind=mm_wp), INTENT(in)           :: p_prod
      !!  Aerosol production pressure level in Pa.
    REAL(kind=mm_wp), INTENT(in)           :: tx_prod
      !! Spherical aerosol mode production rate in \(kg.m^{-2}.s^{-1}\).
    REAL(kind=mm_wp), INTENT(in)           :: rc_prod
      !! Spherical mode characteristic radius for production in meter.
    REAL(kind=mm_wp), INTENT(in)           :: rplanet
      !! Planet radius in meter
    REAL(kind=mm_wp), INTENT(in)           :: g0
      !! Planet gravity acceleration at ground level in \(m.s^{-2}\).
    REAL(kind=mm_wp), INTENT(in)           :: air_rad
      !! Air molecules mean radius in meter.
    REAL(kind=mm_wp), INTENT(in)           :: air_mmol
      !! Air molecules mean molar mass in \(kg.mol^{-1}\).
    LOGICAL, INTENT(in)                    :: clouds
      !! Clouds microphysics control flag.
    CHARACTER(len=*), INTENT(in), OPTIONAL :: cfgpath
      !! Internal microphysic configuration file.

    INTEGER                                           :: coag_choice
    REAL(kind=mm_wp)                                  :: fiad_max,fiad_min,df,rm,rho_aer
    LOGICAL                                           :: w_h_prod,w_h_sed,w_h_coag,w_c_sed,w_c_nucond, &
                                                         no_fiadero,fwsed_m0,fwsed_m3
    TYPE(error)                                       :: err
    INTEGER                                           :: i
    TYPE(cfgparser)                                   :: cparser
    CHARACTER(len=st_slen)                            :: spcpath,pssfile,mqfile,opt_file
    CHARACTER(len=st_slen), DIMENSION(:), ALLOCATABLE :: species
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE       :: tmp
    REAL(kind=mm_wp)                                  :: m0as_min,rcs_min,m0af_min,rcf_min,m0n_min
    LOGICAL                                           :: wdebug

    w_h_prod    = .true.
    w_h_sed     = .true.
    w_h_coag    = .true.
    w_c_sed     = clouds
    w_c_nucond  = clouds
    fwsed_m0    = .true.
    fwsed_m3    = .false.
    no_fiadero  = .true.
    fiad_min    = 0.1_mm_wp
    fiad_max    = 10._mm_wp
    coag_choice = 7
    wdebug      = .false.
    m0as_min    = 1e-10_mm_wp
    rcs_min     = 1e-9_mm_wp
    m0af_min    = 1e-10_mm_wp
    rcf_min     = 1e-9_mm_wp
    m0n_min     = 1e-10_mm_wp

    WRITE(*,'(a)') "##### MMP_GCM SPEAKING #####"
    WRITE(*,'(a)') "I will initialize the microphysics model in moments YAMMS"
    WRITE(*,'(a)') "On error I will simply abort the program. Stay near your computer !"
    WRITE(*,*)
    WRITE(*,'(a)') "Reading muphys configuration file ("//trim(cfgpath)//")..."
    err = cfg_read_config(cparser,TRIM(cfgpath),.true.)
    IF (err /= 0) THEN
      ! RETURN AN ERROR !!
      call abort_program(err)
    ENDIF

    ! YAMMS internal parameters:
    err = mm_check_opt(cfg_get_value(cparser,"rm",rm),rm,50e-9_mm_wp,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"df",df),df,2._mm_wp,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"rho_aer",rho_aer),rho_aer,1000._mm_wp,mm_log)
    ! the following parameters are primarily used to test and debug YAMMS.
    ! They are set in an optional configuration file and default to suitable values for production runs.
    err = mm_check_opt(cfg_get_value(cparser,"haze_production",w_h_prod)          ,w_h_prod   ,.true.     ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"haze_sedimentation",w_h_sed)        ,w_h_sed    ,.true.     ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"haze_coagulation",w_h_coag)         ,w_h_coag   ,.true.     ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"clouds_sedimentation",w_c_sed)      ,w_c_sed    ,clouds     ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"clouds_nucl_cond",w_c_nucond)       ,w_c_nucond ,clouds     ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"wsed_m0",fwsed_m0)                  ,fwsed_m0   ,.true.     ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"wsed_m3",fwsed_m3)                  ,fwsed_m3   ,.false.    ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"no_fiadero",no_fiadero)             ,no_fiadero ,.true.     ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"fiadero_min_ratio",fiad_min)        ,fiad_min   ,0.1_mm_wp  ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"fiadero_max_ratio",fiad_max)        ,fiad_max   ,10._mm_wp  ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"haze_coag_interactions",coag_choice),coag_choice,7          ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"m0as_min",m0as_min)                 ,m0as_min   ,1e-10_mm_wp,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"rcs_min",rcs_min)                   ,rcs_min    ,1e-9_mm_wp ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"m0af_min",m0af_min)                 ,m0af_min   ,1e-10_mm_wp,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"rcf_min",rcf_min)                   ,rcf_min    ,rm         ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"m0n_min",m0n_min)                   ,m0n_min    ,1e-10_mm_wp,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"debug",wdebug)                      ,wdebug     ,.false.    ,mm_log)

    ! Retrieve clouds species configuration file
    spcpath = ''
    IF (clouds) THEN
      err = mm_check_opt(cfg_get_value(cparser,"specie_cfg",spcpath), spcpath, wlog=mm_log)
      IF (err/=0) call abort_program(err)
    ENDIF

    ! Setup alpha function: THEY ARE REQUIRED IN YAMMS global initialization !
    ! spherical mode inter-moments function parameters
    IF (.NOT.cfg_has_section(cparser,'alpha_s')) call abort_program(error("Cannot find [alpha_s] section",-1))
    err = read_aprm(cparser,'alpha_s',mmp_asp)
    IF (err /= 0) call abort_program(error("alpha_s: "//TRIM(err%msg),-1))
    ! fractal mode inter-moments function parameters
    IF (.NOT.cfg_has_section(cparser,'alpha_f')) call abort_program(error("Cannot find [alpha_f] section",-1))
    err = read_aprm(cparser,'alpha_f',mmp_afp)
    IF (err /= 0) call abort_program(error("alpha_s: "//TRIM(err%msg),-1))

    ! YAMMS initialization.
    err = mm_global_init_0(dt,df,rm,rho_aer,p_prod,tx_prod,rc_prod,rplanet,g0, &
                           air_rad,air_mmol,coag_choice,clouds,spcpath,  &
                           w_h_prod,w_h_sed,w_h_coag,w_c_nucond,  &
                           w_c_sed,fwsed_m0,fwsed_m3, &
                           no_fiadero,fiad_min,fiad_max, &
                           m0as_min,rcs_min,m0af_min,rcf_min,m0n_min,wdebug)
    IF (err /= 0) call abort_program(err)

    ! Extra initialization (needed for YAMMS method interfaces)
    err = mm_check_opt(cfg_get_value(cparser, "transfert_probability", mmp_w_ps2s), mmp_w_ps2s, wlog=mm_log)
    IF (err/=0) call abort_program(err)
    err = mm_check_opt(cfg_get_value(cparser, "electric_charging"    , mmp_w_qe  ), mmp_w_qe, wlog=mm_log)
    IF (err/=0) call abort_program(err)

    ! initialize transfert probabilities look-up tables
    IF (mm_w_haze_coag .AND. mmp_w_ps2s) THEN
      err = mm_check_opt(cfg_get_value(cparser, "ps2s_file", pssfile), pssfile)
      IF (err /= 0) call abort_program(err)

      IF (.NOT.read_dset(pssfile,'p_m0_co',mmp_pco0p)) THEN
        call abort_program(error("Cannot get 'p_m0_co' from "//pssfile,-1))
      ENDIF
      IF (.NOT.read_dset(pssfile,'p_m3_co',mmp_pco3p)) THEN
        call abort_program(error("Cannot get 'p_m3_co' from "//pssfile,-1))
      ENDIF
      IF (.NOT.read_dset(pssfile,'p_m0_fm',mmp_pfm0p)) THEN
        call abort_program(error("Cannot get 'p_m0_fm' from "//pssfile,-1))
      ENDIF
      IF (.NOT.read_dset(pssfile,'p_m3_fm',mmp_pfm3p)) THEN
        call abort_program(error("Cannot get 'p_m3_fm' from "//pssfile,-1))
      ENDIF
    ENDIF
    ! initialize mean electric correction look-up tables
    IF (mm_w_haze_coag .AND. mmp_w_qe) THEN
      err = mm_check_opt(cfg_get_value(cparser, "mq_file", mqfile), mqfile)
      IF (err /= 0) call abort_program(err)

      IF (.NOT.read_dset(mqfile,'qbsf0',mmp_qbsf0)) THEN
        call abort_program(error("Cannot get 'qbsf0' from "//mqfile,-1))
      ELSE
        mmp_qbsf0_e(1,1) = MINVAL(mmp_qbsf0%x)
        mmp_qbsf0_e(1,2) = MAXVAL(mmp_qbsf0%x)
        mmp_qbsf0_e(2,1) = MINVAL(mmp_qbsf0%y)
        mmp_qbsf0_e(2,2) = MAXVAL(mmp_qbsf0%y)
      ENDIF
      IF (.NOT.read_dset(mqfile,'qbsf3',mmp_qbsf3)) THEN
        call abort_program(error("Cannot get 'qbsf3' from "//mqfile,-1))
      ELSE
        mmp_qbsf3_e(1,1) = MINVAL(mmp_qbsf3%x)
        mmp_qbsf3_e(1,2) = MAXVAL(mmp_qbsf3%x)
        mmp_qbsf3_e(2,1) = MINVAL(mmp_qbsf3%y)
        mmp_qbsf3_e(2,2) = MAXVAL(mmp_qbsf3%y)
      ENDIF
      IF (.NOT.read_dset(mqfile,'qbff0',mmp_qbff0)) THEN
        call abort_program(error("Cannot get 'qbff0' from "//mqfile,-1))
      ELSE
        mmp_qbff0_e(1,1) = MINVAL(mmp_qbff0%x)
        mmp_qbff0_e(1,2) = MAXVAL(mmp_qbff0%x)
        mmp_qbff0_e(2,1) = MINVAL(mmp_qbff0%y)
        mmp_qbff0_e(2,2) = MAXVAL(mmp_qbff0%y)
      ENDIF
    ENDIF

    ! btk coefficients
    IF (.NOT.cfg_has_section(cparser,'btks')) call abort_program(error("Cannot find [btks] section",-1))
    err = cfg_get_value(cparser,"btks/bt0",tmp) ; IF (err/=0) call abort_program(error("bt0: "//TRIM(err%msg),-1))
    IF (SIZE(tmp) /= 5)  call abort_program(error("bt0: Inconsistent number of coefficients",-1))
    mmp_bt0 = tmp
    err = cfg_get_value(cparser,"btks/bt3",tmp) ; IF (err/=0) call abort_program(error("bt3: "//TRIM(err%msg),-1))
    IF (SIZE(tmp) /= 5)  call abort_program(error("bt3: Inconsistent number of coefficients",-1))
    mmp_bt3 = tmp

    ! dump parameters ...
    WRITE(*,'(a)')        "========= MUPHYS PARAMETERS ==========="
    WRITE(*,'(a,L2)')     "transfert_probability: ", mmp_w_ps2s
    WRITE(*,'(a,L2)')     "electric_charging    : ", mmp_w_qe
    call mm_dump_parameters()
    IF (clouds) THEN
      DO i=1, size(mm_xESPS)
        print*, TRIM(mm_xESPS(i)%name), " fmol2fmas = ", mm_xESPS(i)%fmol2fmas
      ENDDO
    ENDIF

  END SUBROUTINE mmp_initialize

  FUNCTION read_aprm(parser,sec,pp) RESULT(err)
    !! Read and store [[mmp_gcm(module):aprm(type)]] parameters.
    TYPE(cfgparser), INTENT(in)  :: parser !! Configuration parser
    CHARACTER(len=*), INTENT(in) :: sec    !! Name of the section that contains the parameters.
    TYPE(aprm), INTENT(out)      :: pp     !! [[mmp_gcm(module):aprm(type)]] object that stores the parameters values.
    TYPE(error) :: err                     !! Error status of the function.
    err = cfg_get_value(parser,TRIM(sec)//'/a',pp%a) ; IF (err /= 0) RETURN
    err = cfg_get_value(parser,TRIM(sec)//'/b',pp%b) ; IF (err /= 0) RETURN
    err = cfg_get_value(parser,TRIM(sec)//'/c',pp%c) ; IF (err /= 0) RETURN
    IF (SIZE(pp%a) /= SIZE(pp%b) .OR. SIZE(pp%a) /= SIZE(pp%c)) &
    err = error("Inconsistent number of coefficients (a,b, and c must have the same size)",-1)
    RETURN
  END FUNCTION read_aprm

END MODULE MMP_GCM