Ignore:
Timestamp:
Oct 12, 2023, 10:30:22 AM (15 months ago)
Author:
slebonnois
Message:

BBT : Update for the titan microphysics and cloud model

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/libf/muphytitan/mmp_gcm.f90

    r1926 r3083  
    1 ! Copyright 2017 Université de Reims Champagne-Ardenne 
     1! Copyright 2017 Université de Reims Champagne-Ardenne
    22! Contributor: J. Burgalat (GSMA, URCA)
    33! email of the author : jeremie.burgalat@univ-reims.fr
    4 ! 
     4!
    55! This software is a computer program whose purpose is to compute
    66! microphysics processes using a two-moments scheme.
    7 ! 
     7!
    88! This library is governed by the CeCILL license under French law and
    9 ! abiding by the rules of distribution of free software.  You can  use, 
     9! abiding by the rules of distribution of free software.  You can  use,
    1010! modify and/ or redistribute the software under the terms of the CeCILL
    1111! license as circulated by CEA, CNRS and INRIA at the following URL
    12 ! "http://www.cecill.info". 
    13 ! 
     12! "http://www.cecill.info".
     13!
    1414! As a counterpart to the access to the source code and  rights to copy,
    1515! modify and redistribute granted by the license, users are provided only
    1616! with a limited warranty  and the software's author,  the holder of the
    1717! economic rights,  and the successive licensors  have only  limited
    18 ! liability. 
    19 ! 
     18! liability.
     19!
    2020! In this respect, the user's attention is drawn to the risks associated
    2121! with loading,  using,  modifying and/or developing or reproducing the
    … …  
    2525! professionals having in-depth computer knowledge. Users are therefore
    2626! encouraged to load and test the software's suitability as regards their
    27 ! requirements in conditions enabling the security of their systems and/or 
    28 ! data to be ensured and,  more generally, to use and operate it in the 
    29 ! same conditions as regards security. 
    30 ! 
     27! requirements in conditions enabling the security of their systems and/or
     28! data to be ensured and,  more generally, to use and operate it in the
     29! same conditions as regards security.
     30!
    3131! The fact that you are presently reading this means that you have had
    3232! knowledge of the CeCILL license and that you accept its terms.
    … …  
    4444  IMPLICIT NONE
    4545
    46   CONTAINS 
    47    
     46  CONTAINS
     47
    4848  SUBROUTINE mmp_initialize(dt,p_prod,tx_prod,rc_prod,rplanet,g0, air_rad,air_mmol,clouds,cfgpath)
    4949    !! Initialize global parameters of the model.
    50     !! 
     50    !!
    5151    !! The function initializes all the global parameters of the model from direct input.
    52     !! Boolean (and Fiadero) parameters are optional as they are rather testing parameters. Their 
    53     !! default values are suitable for production runs. 
     52    !! Boolean (and Fiadero) parameters are optional as they are rather testing parameters. Their
     53    !! default values are suitable for production runs.
    5454    !! @note
    5555    !! If the subroutine fails to initialize parameters, the run is aborted.
    5656    !!
    5757    !! @warning
    58     !! If OpenMP is activated, this subroutine must be called in an $OMP SINGLE statement as it 
     58    !! If OpenMP is activated, this subroutine must be called in an $OMP SINGLE statement as it
    5959    !! initializes global variable that are not thread private.
    6060    !!
    … …  
    8787
    8888    INTEGER                                           :: coag_choice
    89     REAL(kind=mm_wp)                                  :: fiad_max, fiad_min,df,rm,rho_aer
    90     LOGICAL                                           :: w_h_prod, w_h_sed, w_h_coag, w_c_sed, w_c_nucond, &
    91                                                          no_fiadero, fwsed_m0, fwsed_m3
     89    REAL(kind=mm_wp)                                  :: fiad_max,fiad_min,df,rm,rho_aer
     90    LOGICAL                                           :: w_h_prod,w_h_sed,w_h_coag,w_c_sed,w_c_nucond, &
     91                                                         no_fiadero,fwsed_m0,fwsed_m3
    9292    TYPE(error)                                       :: err
    9393    INTEGER                                           :: i
    … …  
    9696    CHARACTER(len=st_slen), DIMENSION(:), ALLOCATABLE :: species
    9797    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE       :: tmp
     98    REAL(kind=mm_wp)                                  :: m0as_min,rcs_min,m0af_min,rcf_min,m0n_min
     99    LOGICAL                                           :: wdebug
    98100
    99101    w_h_prod    = .true.
    … …  
    104106    fwsed_m0    = .true.
    105107    fwsed_m3    = .false.
    106     no_fiadero  = .false.
     108    no_fiadero  = .true.
    107109    fiad_min    = 0.1_mm_wp
    108110    fiad_max    = 10._mm_wp
    109111    coag_choice = 7
     112    wdebug      = .false.
     113    m0as_min    = 1e-10_mm_wp
     114    rcs_min     = 1e-9_mm_wp
     115    m0af_min    = 1e-10_mm_wp
     116    rcf_min     = 1e-9_mm_wp
     117    m0n_min     = 1e-10_mm_wp
    110118
    111119    WRITE(*,'(a)') "##### MMP_GCM SPEAKING #####"
    112     WRITE(*,'(a)') "I will initialize ze microphysics model in moments YAMMS"
     120    WRITE(*,'(a)') "I will initialize the microphysics model in moments YAMMS"
    113121    WRITE(*,'(a)') "On error I will simply abort the program. Stay near your computer !"
    114122    WRITE(*,*)
    115123    WRITE(*,'(a)') "Reading muphys configuration file ("//trim(cfgpath)//")..."
    116     err = cfg_read_config(cparser,TRIM(cfgpath),.true.) 
     124    err = cfg_read_config(cparser,TRIM(cfgpath),.true.)
    117125    IF (err /= 0) THEN
    118126      ! RETURN AN ERROR !!
    119127      call abort_program(err)
    120128    ENDIF
    121    
     129
    122130    ! YAMMS internal parameters:
    123131    err = mm_check_opt(cfg_get_value(cparser,"rm",rm),rm,50e-9_mm_wp,mm_log)
    … …  
    126134    ! the following parameters are primarily used to test and debug YAMMS.
    127135    ! They are set in an optional configuration file and default to suitable values for production runs.
    128     err = mm_check_opt(cfg_get_value(cparser,"haze_production",w_h_prod)    ,w_h_prod   ,.true.   ,mm_log)
    129     err = mm_check_opt(cfg_get_value(cparser,"haze_sedimentation",w_h_sed)  ,w_h_sed    ,.true.   ,mm_log)
    130     err = mm_check_opt(cfg_get_value(cparser,"haze_coagulation",w_h_coag)   ,w_h_coag   ,.true.   ,mm_log)
    131     err = mm_check_opt(cfg_get_value(cparser,"clouds_sedimentation",w_c_sed),w_c_sed    ,clouds   ,mm_log)
    132     err = mm_check_opt(cfg_get_value(cparser,"clouds_nucl_cond",w_c_nucond) ,w_c_nucond ,clouds   ,mm_log)
    133     err = mm_check_opt(cfg_get_value(cparser,"wsed_m0",fwsed_m0)            ,fwsed_m0   ,.true.   ,mm_log)
    134     err = mm_check_opt(cfg_get_value(cparser,"wsed_m3",fwsed_m3)            ,fwsed_m3   ,.false.  ,mm_log)
    135     err = mm_check_opt(cfg_get_value(cparser,"no_fiadero",no_fiadero)       ,no_fiadero ,.false.  ,mm_log)
    136     err = mm_check_opt(cfg_get_value(cparser,"fiadero_min_ratio",fiad_min)  ,fiad_min   ,0.1_mm_wp,mm_log)
    137     err = mm_check_opt(cfg_get_value(cparser,"fiadero_max_ratio",fiad_max)  ,fiad_max   ,10._mm_wp,mm_log)
    138     err = mm_check_opt(cfg_get_value(cparser,"haze_coag_interactions",coag_choice),coag_choice,7,mm_log)
    139 
    140     ! optic look-up table file path.
    141     mmp_optic_file = ''
    142     opt_file = ''
    143     err = mm_check_opt(cfg_get_value(cparser,"optics/optic_file",opt_file),opt_file,'',mm_log)
    144     IF (err /= 0) THEN
    145       WRITE(*,'(a)') "Warning: I was unable to retrieve the path of the optic look-up table file:"
    146       WRITE(*,'(a)') "  The GCM may abort if it uses YAMMS optical properties calculation module !"
    147     ELSE
    148       mmp_optic_file = TRIM(opt_file)
    149     ENDIF
     136    err = mm_check_opt(cfg_get_value(cparser,"haze_production",w_h_prod)          ,w_h_prod   ,.true.     ,mm_log)
     137    err = mm_check_opt(cfg_get_value(cparser,"haze_sedimentation",w_h_sed)        ,w_h_sed    ,.true.     ,mm_log)
     138    err = mm_check_opt(cfg_get_value(cparser,"haze_coagulation",w_h_coag)         ,w_h_coag   ,.true.     ,mm_log)
     139    err = mm_check_opt(cfg_get_value(cparser,"clouds_sedimentation",w_c_sed)      ,w_c_sed    ,clouds     ,mm_log)
     140    err = mm_check_opt(cfg_get_value(cparser,"clouds_nucl_cond",w_c_nucond)       ,w_c_nucond ,clouds     ,mm_log)
     141    err = mm_check_opt(cfg_get_value(cparser,"wsed_m0",fwsed_m0)                  ,fwsed_m0   ,.true.     ,mm_log)
     142    err = mm_check_opt(cfg_get_value(cparser,"wsed_m3",fwsed_m3)                  ,fwsed_m3   ,.false.    ,mm_log)
     143    err = mm_check_opt(cfg_get_value(cparser,"no_fiadero",no_fiadero)             ,no_fiadero ,.true.     ,mm_log)
     144    err = mm_check_opt(cfg_get_value(cparser,"fiadero_min_ratio",fiad_min)        ,fiad_min   ,0.1_mm_wp  ,mm_log)
     145    err = mm_check_opt(cfg_get_value(cparser,"fiadero_max_ratio",fiad_max)        ,fiad_max   ,10._mm_wp  ,mm_log)
     146    err = mm_check_opt(cfg_get_value(cparser,"haze_coag_interactions",coag_choice),coag_choice,7          ,mm_log)
     147    err = mm_check_opt(cfg_get_value(cparser,"m0as_min",m0as_min)                 ,m0as_min   ,1e-10_mm_wp,mm_log)
     148    err = mm_check_opt(cfg_get_value(cparser,"rcs_min",rcs_min)                   ,rcs_min    ,1e-9_mm_wp ,mm_log)
     149    err = mm_check_opt(cfg_get_value(cparser,"m0af_min",m0af_min)                 ,m0af_min   ,1e-10_mm_wp,mm_log)
     150    err = mm_check_opt(cfg_get_value(cparser,"rcf_min",rcf_min)                   ,rcf_min    ,rm         ,mm_log)
     151    err = mm_check_opt(cfg_get_value(cparser,"m0n_min",m0n_min)                   ,m0n_min    ,1e-10_mm_wp,mm_log)
     152    err = mm_check_opt(cfg_get_value(cparser,"debug",wdebug)                      ,wdebug     ,.false.    ,mm_log)
    150153
    151154    ! Retrieve clouds species configuration file
    … …  
    156159    ENDIF
    157160
    158     ! YAMMS initialization.
    159     err = mm_global_init_0(dt,df,rm,rho_aer,p_prod,tx_prod,rc_prod,rplanet,g0, &
    160                            air_rad,air_mmol,coag_choice,clouds,spcpath,  &
    161                            w_h_prod,w_h_sed,w_h_coag,w_c_nucond,  &
    162                            w_c_sed,fwsed_m0,fwsed_m3, &
    163                            no_fiadero,fiad_min,fiad_max)
    164     IF (err /= 0) call abort_program(err)
    165 
    166     ! Extra initialization (needed for YAMMS method interfaces)
    167     err = mm_check_opt(cfg_get_value(cparser, "transfert_probability", mmp_w_ps2s), mmp_w_ps2s, wlog=mm_log)
    168     IF (err/=0) call abort_program(err)
    169     err = mm_check_opt(cfg_get_value(cparser, "electric_charging"    , mmp_w_qe  ), mmp_w_qe, wlog=mm_log)
    170     IF (err/=0) call abort_program(err)
    171 
    172     ! initialize transfert probabilities look-up tables
    173     IF (mm_w_haze_coag .AND. mmp_w_ps2s) THEN
    174       err = mm_check_opt(cfg_get_value(cparser, "ps2s_file", pssfile), pssfile)
    175       IF (err /= 0) call abort_program(err)
    176 
    177       IF (.NOT.read_dset(pssfile,'p_m0_co',mmp_pco0p)) THEN
    178         call abort_program(error("Cannot get 'p_m0_co' from "//pssfile,-1))
    179       ENDIF
    180       IF (.NOT.read_dset(pssfile,'p_m3_co',mmp_pco3p)) THEN
    181         call abort_program(error("Cannot get 'p_m3_co' from "//pssfile,-1))
    182       ENDIF
    183       IF (.NOT.read_dset(pssfile,'p_m0_fm',mmp_pfm0p)) THEN
    184         call abort_program(error("Cannot get 'p_m0_fm' from "//pssfile,-1))
    185       ENDIF
    186       IF (.NOT.read_dset(pssfile,'p_m3_fm',mmp_pfm3p)) THEN
    187         call abort_program(error("Cannot get 'p_m3_fm' from "//pssfile,-1))
    188       ENDIF
    189     ENDIF
    190     ! initialize mean electric correction look-up tables
    191     IF (mm_w_haze_coag .AND. mmp_w_qe) THEN
    192       err = mm_check_opt(cfg_get_value(cparser, "mq_file", mqfile), mqfile)
    193       IF (err /= 0) call abort_program(err)
    194 
    195       IF (.NOT.read_dset(mqfile,'qbsf0',mmp_qbsf0)) THEN
    196         call abort_program(error("Cannot get 'qbsf0' from "//mqfile,-1))
    197       ELSE
    198         mmp_qbsf0_e(1,1) = MINVAL(mmp_qbsf0%x)
    199         mmp_qbsf0_e(1,2) = MAXVAL(mmp_qbsf0%x)
    200         mmp_qbsf0_e(2,1) = MINVAL(mmp_qbsf0%y)
    201         mmp_qbsf0_e(2,2) = MAXVAL(mmp_qbsf0%y)
    202       ENDIF
    203       IF (.NOT.read_dset(mqfile,'qbsf3',mmp_qbsf3)) THEN
    204         call abort_program(error("Cannot get 'qbsf3' from "//mqfile,-1))
    205       ELSE
    206         mmp_qbsf3_e(1,1) = MINVAL(mmp_qbsf3%x)
    207         mmp_qbsf3_e(1,2) = MAXVAL(mmp_qbsf3%x)
    208         mmp_qbsf3_e(2,1) = MINVAL(mmp_qbsf3%y)
    209         mmp_qbsf3_e(2,2) = MAXVAL(mmp_qbsf3%y)
    210       ENDIF
    211       IF (.NOT.read_dset(mqfile,'qbff0',mmp_qbff0)) THEN
    212         call abort_program(error("Cannot get 'qbff0' from "//mqfile,-1))
    213       ELSE
    214         mmp_qbff0_e(1,1) = MINVAL(mmp_qbff0%x)
    215         mmp_qbff0_e(1,2) = MAXVAL(mmp_qbff0%x)
    216         mmp_qbff0_e(2,1) = MINVAL(mmp_qbff0%y)
    217         mmp_qbff0_e(2,2) = MAXVAL(mmp_qbff0%y)
    218       ENDIF
    219     ENDIF
     161    ! Setup alpha function: THEY ARE REQUIRED IN YAMMS global initialization !
    220162    ! spherical mode inter-moments function parameters
    221163    IF (.NOT.cfg_has_section(cparser,'alpha_s')) call abort_program(error("Cannot find [alpha_s] section",-1))
    … …  
    227169    IF (err /= 0) call abort_program(error("alpha_s: "//TRIM(err%msg),-1))
    228170
    229     ! get size-distribution laws parameters
    230     IF (.NOT.cfg_has_section(cparser,'dndr_s')) call abort_program(error("Cannot find [dndr_s] section",-2))
    231     err = read_nprm(cparser,'dndr_s',mmp_pns)
    232     IF (err /= 0) call abort_program(error("dndr_s: "//TRIM(err%msg),-2))
    233     IF (.NOT.cfg_has_section(cparser,'dndr_f')) call abort_program(error("Cannot find [dndr_f] section",-2))
    234     err = read_nprm(cparser,'dndr_f',mmp_pnf)
    235     IF (err /= 0) call abort_program(error("dndr_f: "//TRIM(err%msg),-2))
     171    ! YAMMS initialization.
     172    err = mm_global_init_0(dt,df,rm,rho_aer,p_prod,tx_prod,rc_prod,rplanet,g0, &
     173                           air_rad,air_mmol,coag_choice,clouds,spcpath,  &
     174                           w_h_prod,w_h_sed,w_h_coag,w_c_nucond,  &
     175                           w_c_sed,fwsed_m0,fwsed_m3, &
     176                           no_fiadero,fiad_min,fiad_max, &
     177                           m0as_min,rcs_min,m0af_min,rcf_min,m0n_min,wdebug)
     178    IF (err /= 0) call abort_program(err)
     179
     180    ! Extra initialization (needed for YAMMS method interfaces)
     181    err = mm_check_opt(cfg_get_value(cparser, "transfert_probability", mmp_w_ps2s), mmp_w_ps2s, wlog=mm_log)
     182    IF (err/=0) call abort_program(err)
     183    err = mm_check_opt(cfg_get_value(cparser, "electric_charging"    , mmp_w_qe  ), mmp_w_qe, wlog=mm_log)
     184    IF (err/=0) call abort_program(err)
     185
     186    ! initialize transfert probabilities look-up tables
     187    IF (mm_w_haze_coag .AND. mmp_w_ps2s) THEN
     188      err = mm_check_opt(cfg_get_value(cparser, "ps2s_file", pssfile), pssfile)
     189      IF (err /= 0) call abort_program(err)
     190
     191      IF (.NOT.read_dset(pssfile,'p_m0_co',mmp_pco0p)) THEN
     192        call abort_program(error("Cannot get 'p_m0_co' from "//pssfile,-1))
     193      ENDIF
     194      IF (.NOT.read_dset(pssfile,'p_m3_co',mmp_pco3p)) THEN
     195        call abort_program(error("Cannot get 'p_m3_co' from "//pssfile,-1))
     196      ENDIF
     197      IF (.NOT.read_dset(pssfile,'p_m0_fm',mmp_pfm0p)) THEN
     198        call abort_program(error("Cannot get 'p_m0_fm' from "//pssfile,-1))
     199      ENDIF
     200      IF (.NOT.read_dset(pssfile,'p_m3_fm',mmp_pfm3p)) THEN
     201        call abort_program(error("Cannot get 'p_m3_fm' from "//pssfile,-1))
     202      ENDIF
     203    ENDIF
     204    ! initialize mean electric correction look-up tables
     205    IF (mm_w_haze_coag .AND. mmp_w_qe) THEN
     206      err = mm_check_opt(cfg_get_value(cparser, "mq_file", mqfile), mqfile)
     207      IF (err /= 0) call abort_program(err)
     208
     209      IF (.NOT.read_dset(mqfile,'qbsf0',mmp_qbsf0)) THEN
     210        call abort_program(error("Cannot get 'qbsf0' from "//mqfile,-1))
     211      ELSE
     212        mmp_qbsf0_e(1,1) = MINVAL(mmp_qbsf0%x)
     213        mmp_qbsf0_e(1,2) = MAXVAL(mmp_qbsf0%x)
     214        mmp_qbsf0_e(2,1) = MINVAL(mmp_qbsf0%y)
     215        mmp_qbsf0_e(2,2) = MAXVAL(mmp_qbsf0%y)
     216      ENDIF
     217      IF (.NOT.read_dset(mqfile,'qbsf3',mmp_qbsf3)) THEN
     218        call abort_program(error("Cannot get 'qbsf3' from "//mqfile,-1))
     219      ELSE
     220        mmp_qbsf3_e(1,1) = MINVAL(mmp_qbsf3%x)
     221        mmp_qbsf3_e(1,2) = MAXVAL(mmp_qbsf3%x)
     222        mmp_qbsf3_e(2,1) = MINVAL(mmp_qbsf3%y)
     223        mmp_qbsf3_e(2,2) = MAXVAL(mmp_qbsf3%y)
     224      ENDIF
     225      IF (.NOT.read_dset(mqfile,'qbff0',mmp_qbff0)) THEN
     226        call abort_program(error("Cannot get 'qbff0' from "//mqfile,-1))
     227      ELSE
     228        mmp_qbff0_e(1,1) = MINVAL(mmp_qbff0%x)
     229        mmp_qbff0_e(1,2) = MAXVAL(mmp_qbff0%x)
     230        mmp_qbff0_e(2,1) = MINVAL(mmp_qbff0%y)
     231        mmp_qbff0_e(2,2) = MAXVAL(mmp_qbff0%y)
     232      ENDIF
     233    ENDIF
    236234
    237235    ! btk coefficients
    … …  
    249247    WRITE(*,'(a,L2)')     "electric_charging    : ", mmp_w_qe
    250248    call mm_dump_parameters()
    251      
     249
     250    IF (clouds) THEN
     251      DO i=1, size(mm_xESPS)
     252        print*, TRIM(mm_xESPS(i)%name), " fmol2fmas = ", mm_xESPS(i)%fmol2fmas
     253      ENDDO
     254    ENDIF
     255
    252256  END SUBROUTINE mmp_initialize
    253257
    254258  FUNCTION read_aprm(parser,sec,pp) RESULT(err)
    255     !! Read and store [[mmp_gcm(module):aprm(type)]] parameters. 
    256     TYPE(cfgparser), INTENT(in)  :: parser !! Configuration parser 
     259    !! Read and store [[mmp_gcm(module):aprm(type)]] parameters.
     260    TYPE(cfgparser), INTENT(in)  :: parser !! Configuration parser
    257261    CHARACTER(len=*), INTENT(in) :: sec    !! Name of the section that contains the parameters.
    258262    TYPE(aprm), INTENT(out)      :: pp     !! [[mmp_gcm(module):aprm(type)]] object that stores the parameters values.
    … …  
    266270  END FUNCTION read_aprm
    267271
    268   FUNCTION read_nprm(parser,sec,pp) RESULT(err)
    269     !! Read and store [[mmp_gcm(module):nprm(type)]] parameters.
    270     TYPE(cfgparser), INTENT(in)  :: parser !! Configuration parser
    271     CHARACTER(len=*), INTENT(in) :: sec    !! Name of the section that contains the parameters.
    272     TYPE(nprm), INTENT(out)      :: pp     !! [[mmp_gcm(module):nprm(type)]] object that stores the parameters values.
    273     TYPE(error) :: err                     !! Error status of the function.
    274     err = cfg_get_value(parser,TRIM(sec)//'/rc',pp%rc) ; IF (err /= 0) RETURN
    275     err = cfg_get_value(parser,TRIM(sec)//'/a0',pp%a0) ; IF (err /= 0) RETURN
    276     err = cfg_get_value(parser,TRIM(sec)//'/c',pp%c)   ; IF (err /= 0) RETURN
    277     err = cfg_get_value(parser,TRIM(sec)//'/a',pp%a)   ; IF (err /= 0) RETURN
    278     err = cfg_get_value(parser,TRIM(sec)//'/b',pp%b)   ; IF (err /= 0) RETURN
    279     IF (SIZE(pp%a) /= SIZE(pp%b)) &
    280       err = error("Inconsistent number of coefficients (a and b must have the same size)",3)
    281     RETURN
    282   END FUNCTION read_nprm
    283 
    284272END MODULE MMP_GCM
    285273
Note: See TracChangeset for help on using the changeset viewer.