Changeset 2242 for trunk/LMDZ.TITAN/libf


Ignore:
Timestamp:
Feb 17, 2020, 5:44:12 PM (6 years ago)
Author:
jvatant
Message:

+ Update the interface with YAMMS so it now correctly handles the small values of the moments,
requiring dynamics to have a threshold quite low (set to 1D-200) and a sanity check in calmufi
corresponding to this value. Thus it removes 'most' of the unphysical radius obtained in
YAMMS. There are still some, but at least there is no more problem of model stability for the moments
and the code can run.
Still, take care the day you want to calculate opacities from the radii and not the moments.
Although, note that there are some negative values, in the output files, but theses negatives are

harmless, as they are only present in output files, dynamics reseting to epsilon after.

+ Given theses pbs with the radii, also update the optics so that it computes the opacities in
a simpler way, directly for M3, through look-up tables, M3 being a good proxy.
--JVO

Location:
trunk/LMDZ.TITAN/libf
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r1897 r2242  
    10661066    ! il faudrait peut etre revoir la gestion des zeros ici et la
    10671067    ! remplacer par une valeur seuil des moments.
     1068    !
     1069    !-> JVO 19 : Done. Zero threshold set at espilon from dynamics on the
     1070    ! input moments in calmufi (safer than here). Might still be some unphysical
     1071    ! values after the dynamics near the threshold. Could be a could idea to add
     1072    ! a sanity check filtering too high values of radii.
     1073    !
     1074    ! TBD : Add a sanity check for high radii ????
    10681075    WHERE(mm_m3aer_s > 0._mm_wp .AND. mm_m0aer_s > 0._mm_wp)
    10691076      mm_rcs = mm_get_rcs(mm_m0aer_s,mm_m3aer_s)

  • trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90

    r2109 r2242  
    9292
    9393  ! Initialization of zdq here since intent=out and no action performed on every tracers
    94   zdq(:,:,:) = 0.0
     94  zdq(:,:,:) = 0.D0
    9595
    9696  ! Initialize tracers updated with former processes from physics
    … …  
    102102    ! We suppose a given order of tracers !
    103103    int2ext(ilon,:) = ( plev(ilon,1:nlay)-plev(ilon,2:nlay+1) ) / g3d(ilon,1:nlay)
    104     m0as(:) = zq(ilon,:,1) * int2ext(ilon,:)
    105     m3as(:) = zq(ilon,:,2) * int2ext(ilon,:)
    106     m0af(:) = zq(ilon,:,3) * int2ext(ilon,:)
    107     m3af(:) = zq(ilon,:,4) * int2ext(ilon,:)
     104   
     105    ! Check because of the threshold of small tracers values in the dynamics
     106    ! WARNING : With this patch it enables to handles the small values required
     107    ! by YAMMS, but still it might still leads to some unphysical values of
     108    ! radii inside YAMMS, harmless for now, but who might be a problem the day
     109    ! you'll want to compute optics from the radii.
     110    WHERE (pq(ilon,:,1) > 2.D-200 .AND. pq(ilon,:,2) > 2.D-200) ! Test on both moments to avoid divergences if one hit the threshold but not the other
     111      m0as(:) = zq(ilon,:,1) * int2ext(ilon,:)                  ! It can still be a pb if both m0 and m3 has been set to epsilon at the beginning of dynamics
     112      m3as(:) = zq(ilon,:,2) * int2ext(ilon,:)                  ! then mixed, even though both are above the threshold, their ratio can be nonsense
     113    ELSEWHERE
     114      m0as(:)=0.D0
     115      m3as(:)=0.D0
     116    ENDWHERE
     117    WHERE (pq(ilon,:,3) > 2.D-200 .AND. pq(ilon,:,4) > 2.D-200)
     118      m0af(:) = zq(ilon,:,3) * int2ext(ilon,:)
     119      m3af(:) = zq(ilon,:,4) * int2ext(ilon,:)
     120    ELSEWHERE
     121      m0af(:)=0.D0
     122      m3af(:)=0.D0
     123    ENDWHERE
    108124   
    109125    if (callclouds) then ! if call clouds
    … …  
    131147    !dm0n(:) = 0._mm_wp ; dm3n(:) = 0._mm_wp ; dm3i(:,:) = 0._mm_wp ; dgazs(:,:) = 0._mm_wp
    132148
    133     dm0as(:) = 0.0 ; dm3as(:) = 0.0 ; dm0af(:) = 0.0 ; dm3af(:) = 0.0
    134     dm0n(:) = 0.0 ; dm3n(:) = 0.0 ; dm3i(:,:) = 0.0 ; dgazs(:,:) = 0.0
     149    dm0as(:) = 0.D0 ; dm3as(:) = 0.D0 ; dm0af(:) = 0.D0 ; dm3af(:) = 0.D0
     150    dm0n(:) = 0.D0 ; dm3n(:) = 0.D0 ; dm3i(:,:) = 0.D0 ; dgazs(:,:) = 0.D0
    135151    ! call microphysics
    136152
    … …  
    165181      enddo
    166182    endif
    167 
    168     ! Sanity check ( way safer to be done here rather than within YAMMS )
    169     WHERE (zq+zdq < 0.0) ; zdq = -zq ; END WHERE
    170 
     183   
    171184  END DO ! loop on ilon
    172185

  • trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90

    r1822 r2242  
    1515     
    1616      ! Default directories for correlated-k data
    17       ! Set in inifis_mod
     17      ! Read in inifis_mod, set in physiq_mod
    1818      character(LEN=100),save :: corrkdir = 'datagcm/corrk_data/50X50'
    1919      character(LEN=100),save :: banddir  = 'datagcm/corrk_data/50X50/50x50'
    … …  
    2525!$OMP THREADPRIVATE(config_mufi)
    2626
     27      ! Default file for coupled microphysics optical properties
     28      ! Set in physiq_mod
     29      character(LEN=100),save :: haze_opt_file ='datagcm/HAZE_OPT_23x23.DAT'
     30!$OMP THREADPRIVATE(haze_opt_file)
     31
    2732      end module datafile_mod
    2833!-----------------------------------------------------------------------

  • trunk/LMDZ.TITAN/libf/phytitan/optci.F90

    r2133 r2242  
    66                         tgasref,pfgasref,wnoi,scalep,indi
    77  use gases_h
     8  use datafile_mod, only: haze_opt_file
    89  use comcstfi_mod, only: r
    910  use callkeys_mod, only: continuum,graybody,corrk_recombin,               &
    1011                          callclouds,callmufi,seashaze,uncoupl_optic_haze
    1112  use tracer_h, only : nmicro,nice
    12   use MMP_OPTICS
    1313
    1414  implicit none
    … …  
    8383  ! variables for k in units m^-1
    8484  real*8 dz(L_LEVELS)
    85   !real*8 rho !! see test below
    8685
    8786  integer igas, jgas, ilay
    … …  
    8988  integer interm
    9089
    91   real*8 m0as,m3as,m0af,m3af
    92   real*8 ext_s,sca_s,ssa_s,asf_s
    93   real*8 ext_f,sca_f,ssa_f,asf_f
     90  ! Variables for haze optics
     91  character(len=200) file_path
     92  logical file_ok
     93  integer dumch
     94  real*8  dumwvl
     95
     96  real*8 m3as,m3af
     97  real*8 dtauaer_s,dtauaer_f
     98  real*8,save :: rhoaer_s(L_NSPECTI),ssa_s(L_NSPECTI),asf_s(L_NSPECTI)
     99  real*8,save :: rhoaer_f(L_NSPECTI),ssa_f(L_NSPECTI),asf_f(L_NSPECTI)
     100!$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,asf_s,asf_f)
     101
    94102  logical,save :: firstcall=.true.
    95   !$OMP THREADPRIVATE(firstcall)
     103!$OMP THREADPRIVATE(firstcall)
     104
    96105
    97106  !! AS: to save time in computing continuum (see bilinearbig)
    … …  
    101110  ENDIF
    102111
    103   ! Some initialisation beacause there's a pb with disr_haze at the limits (nw=1)
     112  ! Some initialisation because there's a pb with disr_haze at the limits (nw=1)
    104113  ! I should check this - For now we set vars to zero : better than nans - JVO 2017
    105  
    106   dhaze_t(:,:) = 0.
    107   ssa_t(:,:) = 0.
    108   asf_t(:,:) = 0.
     114  DHAZE_T(:,:) = 0.0
     115  SSA_T(:,:)   = 0.0
     116  ASF_T(:,:)   = 0.0
     117
     118  ! Load tabulated haze optical properties if needed.
     119  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     120  IF (firstcall .AND. callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
     121     OPEN(12,file=TRIM(haze_opt_file),form='formatted') ! The file has been inquired in physiq_mod firstcall
     122     READ(12,*) ! dummy header
     123     DO NW=1,L_NSPECTI
     124       READ(12,*) dumch, dumwvl, rhoaer_f(nw), ssa_f(nw), asf_f(nw), rhoaer_s(nw), ssa_s(nw), asf_s(nw)
     125     ENDDO
     126     CLOSE(12)
     127  ENDIF
    109128
    110129  !=======================================================================
    … …  
    118137  do K=2,L_LEVELS
    119138 
    120      ilay = k / 2 ! int. arithmetic => gives the gcm layer index
     139     ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
    121140     
    122141     DPR(k) = PLEV(K)-PLEV(K-1)
    … …  
    137156     do K=2,L_LEVELS
    138157     
    139         ilay = k / 2 ! int. arithmetic => gives the gcm layer index
     158        ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
    140159       
    141         ! Optical coupling of YAMMS is plugged but inactivated for now
    142         ! as long as the microphysics only isn't fully debugged -- JVO 01/18
    143160        IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
    144           m0as = pqmo(ilay,1)
    145           m3as = pqmo(ilay,2)
    146           m0af = pqmo(ilay,3)
    147           m3af = pqmo(ilay,4)
    148 
    149           IF (.NOT.mmp_sph_optics_ir(m0as,m3as,nw,ext_s,sca_s,ssa_s,asf_s)) &
    150           CALL abort_gcm("optcv", "Fatal error in mmp_sph_optics_ir", 12)
    151           IF (.NOT.mmp_fra_optics_ir(m0af,m3af,nw,ext_f,sca_f,ssa_f,asf_f)) &
    152           CALL abort_gcm("optcv", "Fatal error in mmp_fra_optics_ir", 12)
    153           dhaze_T(k,nw) = ext_s+ext_f
    154           SSA_T(k,nw)   = (sca_s+sca_f)/dhaze_T(k,nw)
    155           ASF_T(k,nw)   = (asf_s*sca_s + asf_f*sca_f) /(sca_s+sca_f)
     161          m3as = pqmo(ilay,2) / 2.0
     162          m3af = pqmo(ilay,4) / 2.0
     163         
     164          IF ( ilay .lt. 18 ) THEN
     165            m3as = pqmo(18,2) / 2.0 *dz(k)/dz(18)
     166            m3af = pqmo(18,4) / 2.0 *dz(k)/dz(18)
     167          ENDIF
     168
     169          dtauaer_s     = m3as*rhoaer_s(nw)
     170          dtauaer_f     = m3af*rhoaer_f(nw)
     171          DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
     172         
     173          IF ( dtauaer_s + dtauaer_f .GT. 1.D-30 ) THEN
     174            SSA_T(k,nw)   = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
     175            ASF_T(k,nw)   = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) )  &
     176                            / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f )
     177          ELSE
     178             DHAZE_T(k,nw) = 0.D0
     179             SSA_T(k,nw)   = 1.0
     180             ASF_T(k,nw)   = 1.0
     181          ENDIF
     182         
    156183          IF (callclouds.and.firstcall) &
    157184            WRITE(*,*) 'WARNING: In optci, optical properties &
    … …  
    159186        ELSE
    160187          ! Call fixed vertical haze profile of extinction - same for all columns
    161           call disr_haze(dz(k),plev(k),wnoi(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
    162           if (seashaze) dhaze_T(k,nw) = dhaze_T(k,nw)*seashazefact(k)
     188          call disr_haze(dz(k),plev(k),wnoi(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
     189          if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
    163190        ENDIF
    164191

  • trunk/LMDZ.TITAN/libf/phytitan/optcv.F90

    r2133 r2242  
    66                         tgasref,pfgasref,wnov,scalep,indv
    77  use gases_h
     8  use datafile_mod, only: haze_opt_file
    89  use comcstfi_mod, only: r
    910  use callkeys_mod, only: continuum,graybody,callgasvis,corrk_recombin,     &
    … …  
    99100  integer interm
    100101
    101   real*8 m0as,m3as,m0af,m3af
    102   real*8 ext_s,sca_s,ssa_s,asf_s
    103   real*8 ext_f,sca_f,ssa_f,asf_f
     102  ! Variables for haze optics
     103  character(len=200) file_path
     104  logical file_ok
     105  integer dumch
     106  real*8  dumwvl
     107
     108  real*8 m3as,m3af
     109  real*8 dtauaer_s,dtauaer_f
     110  real*8,save :: rhoaer_s(L_NSPECTV),ssa_s(L_NSPECTV),asf_s(L_NSPECTV)
     111  real*8,save :: rhoaer_f(L_NSPECTV),ssa_f(L_NSPECTV),asf_f(L_NSPECTV)
     112!$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,asf_s,asf_f)
     113 
    104114  logical,save :: firstcall=.true.
    105   !$OMP THREADPRIVATE(firstcall)
     115!$OMP THREADPRIVATE(firstcall)
    106116
    107117
    … …  
    114124  ! Some initialisation beacause there's a pb with disr_haze at the limits (nw=1)
    115125  ! I should check this - For now we set vars to zero : better than nans - JVO 2017
    116  
    117   dhaze_t(:,:) = 0.
    118   ssa_t(:,:) = 0.
    119   asf_t(:,:) = 0.
    120 
     126  DHAZE_T(:,:) = 0.0
     127  SSA_T(:,:)   = 0.0
     128  ASF_T(:,:)   = 0.0
     129 
     130  ! Load tabulated haze optical properties if needed.
     131  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     132  IF (firstcall .AND. callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
     133     OPEN(12,file=TRIM(haze_opt_file),form='formatted') ! The file has been inquired in physiq_mod firstcall
     134     READ(12,*) ! dummy header
     135     DO NW=1,L_NSPECTI
     136       READ(12,*) ! there's IR 1st
     137     ENDDO
     138     DO NW=1,L_NSPECTV
     139       READ(12,*) dumch, dumwvl, rhoaer_f(nw), ssa_f(nw), asf_f(nw), rhoaer_s(nw), ssa_s(nw), asf_s(nw)
     140     ENDDO
     141     CLOSE(12)
     142  ENDIF
    121143
    122144  !=======================================================================
    … …  
    132154  do K=2,L_LEVELS
    133155 
    134      ilay = k / 2 ! int. arithmetic => gives the gcm layer index
     156     ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
    135157 
    136158     DPR(k) = PLEV(K)-PLEV(K-1)
    … …  
    150172  ! Rayleigh scattering
    151173  do NW=1,L_NSPECTV
    152      TRAY(1:4,NW)   = 1d-30
     174     TRAY(1:4,NW)   = 1.d-30
    153175     do K=5,L_LEVELS
    154176        TRAY(K,NW)   = TAURAY(NW) * DPR(K)
    … …  
    159181  do K=2,L_LEVELS
    160182 
    161      ilay = k / 2 ! int. arithmetic => gives the gcm layer index
     183     ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
    162184
    163185     do NW=1,L_NSPECTV
    164186     
    165         ! Optical coupling of YAMMS is plugged but inactivated (if false) for now
    166         ! as long as the microphysics only isn't fully debugged -- JVO 01/18
    167187        IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
    168           m0as = pqmo(ilay,1)
    169           m3as = pqmo(ilay,2)
    170           m0af = pqmo(ilay,3)
    171           m3af = pqmo(ilay,4)
    172 
    173           IF (.NOT.mmp_sph_optics_vis(m0as,m3as,nw,ext_s,sca_s,ssa_s,asf_s)) &
    174           CALL abort_gcm("optcv", "Fatal error in mmp_sph_optics_vis", 12)
    175           IF (.NOT.mmp_fra_optics_vis(m0af,m3af,nw,ext_f,sca_f,ssa_f,asf_f)) &
    176           CALL abort_gcm("optcv", "Fatal error in mmp_fra_optics_vis", 12)
    177           dhaze_T(k,nw) = ext_s+ext_f
    178           SSA_T(k,nw)   = (sca_s+sca_f)/dhaze_T(k,nw)
    179           ASF_T(k,nw)   = (asf_s*sca_s + asf_f*sca_f) /(sca_s+sca_f)
     188          m3as = pqmo(ilay,2) / 2.0
     189          m3af = pqmo(ilay,4) / 2.0
     190         
     191          IF ( ilay .lt. 18 ) THEN
     192            m3as = pqmo(18,2) / 2.0 * dz(k) / dz(18)
     193            m3af = pqmo(18,4) / 2.0 * dz(k) / dz(18)
     194          ENDIF
     195
     196          dtauaer_s     = m3as*rhoaer_s(nw)
     197          dtauaer_f     = m3af*rhoaer_f(nw)
     198          DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
     199
     200          IF ( dtauaer_s + dtauaer_f .GT. 1.D-30 ) THEN
     201            SSA_T(k,nw)   = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
     202            ASF_T(k,nw)   = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) )  &
     203                            / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f )
     204          ELSE
     205             DHAZE_T(k,nw) = 0.D0
     206             SSA_T(k,nw)   = 1.0
     207             ASF_T(k,nw)   = 1.0
     208          ENDIF
     209         
    180210          IF (callclouds.and.firstcall) &
    181211            WRITE(*,*) 'WARNING: In optcv, optical properties &
    … …  
    183213        ELSE
    184214          ! Call fixed vertical haze profile of extinction - same for all columns
    185           call disr_haze(dz(k),plev(k),wnov(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
    186           if (seashaze) dhaze_T(k,nw) = dhaze_T(k,nw)*seashazefact(k)
     215          call disr_haze(dz(k),plev(k),wnov(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
     216          if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
    187217        ENDIF
    188218       
    189219        !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
    190         !   but visible does not handle very well diffusion in first layer.
    191         !   The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise)
    192         !   in the 4 first semilayers in optcv, but not optci.
    193         !   This solves random variations of the sw heating at the model top.
    194         if (k<5)  dhaze_T(K,:) = 0.0
     220        !   but visible does not handle very well diffusion in first layer.
     221        !   The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise)
     222        !   in the 4 first semilayers in optcv, but not optci.
     223        !   This solves random variations of the sw heating at the model top.
     224        if (k<5)  DHAZE_T(K,:) = 0.0
    195225         
    196226        DRAYAER = TRAY(K,NW)
    … …  
    316346  ! Haze scattering
    317347            !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
    318             !   but not in the visible
    319             !   The dhaze_s is thus set to 0 in the 4 first semilayers in optcv, but not optci.
    320             !   This solves random variations of the sw heating at the model top.
     348            !   but not in the visible
     349            !   The dhaze_s is thus set to 0 in the 4 first semilayers in optcv, but not optci.
     350            !   This solves random variations of the sw heating at the model top.
    321351  DO NW=1,L_NSPECTV
    322       DHAZES_T(1:4,NW) = 0.d0
     352    DHAZES_T(1:4,NW) = 0.d0
    323353    DO K=5,L_LEVELS
    324354      DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW) ! effect of scattering albedo on haze
    … …  
    386416  END DO                 ! end full gauss loop
    387417
    388 
    389418  if(firstcall) firstcall = .false.
    390419

  • trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90

    r2138 r2242  
    2121      use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
    2222      use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
    23       use datafile_mod, only: datadir, corrkdir, banddir
     23      use datafile_mod, only: datadir, corrkdir, banddir, haze_opt_file
    2424      use geometry_mod, only: latitude, longitude, cell_area
    2525      USE comgeomfi_h, only: totarea, totarea_planet
    … …  
    402402      real,dimension(ngrid,nlayer,nq) :: dtpq ! (temporary in 01/18)
    403403
     404      logical file_ok
    404405
    405406!-----------------------------------------------------------------------------
    … …  
    417418        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: g3d   !! Latitude-Altitude depending gravitational acceleration (m.s-2).
    418419        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp  !! Temperature at the center of each layer (K). 
    419         REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq    !! Tracers (\(kg.kg^{-1}}\)).
    420         REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: zdqfi !! Tendency from former processes for tracers (\(kg.kg^{-1}}\)).
    421         REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq   !! Microphysical tendency for tracers (\(kg.kg^{-1}}\)).
     420        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq    !! Tracers (\(X.kg^{-1}}\)).
     421        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: zdqfi !! Tendency from former processes for tracers (\(X.kg^{-1}}\)).
     422        REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq   !! Microphysical tendency for tracers (\(X.kg^{-1}}\)).
    422423      END SUBROUTINE calmufi
    423424    END INTERFACE
    … …  
    484485!        Variables set to 0
    485486!        ~~~~~~~~~~~~~~~~~~
    486          dtrad(:,:) = 0.0
    487          fluxrad(:) = 0.0
    488          zdtsw(:,:) = 0.0
    489          zdtlw(:,:) = 0.0
     487         dtrad(:,:) = 0.D0
     488         fluxrad(:) = 0.D0
     489         zdtsw(:,:) = 0.D0
     490         zdtlw(:,:) = 0.D0
    490491
    491492!        Initialize setup for correlated-k radiative transfer
    492 !        JVO 17 : Was in callcorrk firstcall, but we need spectral intervals for microphysics
    493 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     493!        JVO 17 : Was in callcorrk firstcall, but we need spectral intervals for microphysics.
     494!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    494495
    495496         if (corrk) then
    … …  
    512513           int_dtauv(:,:,:) = 0.D0
    513514           
     515           IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
     516             haze_opt_file=trim(datadir)//'/HAZE_OPTIC_'//trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))//'.DAT'
     517             inquire(file=trim(haze_opt_file),exist=file_ok)
     518             if(.not.file_ok) then
     519               write(*,*) 'The file ',TRIM(haze_opt_file),' with the haze optical properties'
     520               write(*,*) 'was not found by optci.F90 ! Check in ', TRIM(datadir)
     521               write(*,*) 'that you have the one corresponding to the given spectral resolution !!'
     522               write(*,*) 'Meanwhile I abort ...'
     523               call abort
     524              endif
     525           ENDIF
     526           
    514527         endif
    515528
    … …  
    536549
    537550         IF ( callmufi ) THEN
    538            ! WARNING JB (27/03/2018): inimufi AND mmp_initialize_optics should be wrapped in an OMP SINGLE statement.
     551           ! WARNING JB (27/03/2018): inimufi should be wrapped in an OMP SINGLE statement.
    539552           call inimufi(ptimestep)
    540            ! Optical coupling of YAMMS is plugged but inactivated for now
    541            ! as long as the microphysics only isn't fully debugged -- JVO 01/18
    542            ! NOTE JB (03/18): mmp_optic_file is initialized in inimufi => mmp_initialize_optics must be called AFTER inimufi
    543            IF (.NOT.uncoupl_optic_haze) call mmp_initialize_optics(mmp_optic_file)
    544553
    545554           ! initialize microphysics diagnostics arrays.
    … …  
    547556
    548557         ENDIF
    549 
    550558
    551559#ifdef CPP_XIOS
    … …  
    582590!        Initialize albedo calculation.
    583591!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    584          albedo(:,:)=0.0
    585           albedo_bareground(:)=0.0
     592         albedo(:,:)=0.D0
     593          albedo_bareground(:)=0.D0
    586594         call surfini(ngrid,nq,qsurf,albedo,albedo_bareground)
    587595         
    … …  
    660668      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    661669
    662       pdt(:,:)    = 0.0     
    663       zdtsurf(:)  = 0.0
    664       pdq(:,:,:)  = 0.0
    665       dqsurf(:,:) = 0.0
    666       pdu(:,:)    = 0.0
    667       pdv(:,:)    = 0.0
    668       pdpsrf(:)   = 0.0
    669       zflubid(:)  = 0.0     
    670       taux(:)     = 0.0
    671       tauy(:)     = 0.0
     670      pdt(:,:)    = 0.D0     
     671      zdtsurf(:)  = 0.D0
     672      pdq(:,:,:)  = 0.D0
     673      dqsurf(:,:) = 0.D0
     674      pdu(:,:)    = 0.D0
     675      pdv(:,:)    = 0.D0
     676      pdpsrf(:)   = 0.D0
     677      zflubid(:)  = 0.D0     
     678      taux(:)     = 0.D0
     679      tauy(:)     = 0.D0
    672680
    673681      zday=pday+ptime ! Compute time, in sols (and fraction thereof).
    … …  
    906914               fluxtop_lw(:)     = emis(:)*sigma*tsurf(:)**4
    907915
    908                dtrad(:,:)=0.0 ! no atmospheric radiative heating
     916               dtrad(:,:)=0.D0 ! no atmospheric radiative heating
    909917
    910918            endif ! end of corrk
    … …  
    10441052
    10451053         zdh(:,:) = pdt(:,:)/zpopsk(:,:)
    1046          zduadj(:,:)=0.0
    1047          zdvadj(:,:)=0.0
    1048          zdhadj(:,:)=0.0
     1054         zduadj(:,:)=0.D0
     1055         zdvadj(:,:)=0.D0
     1056         zdhadj(:,:)=0.D0
    10491057
    10501058
    … …  
    10881096 
    10891097         if (callmufi) then
    1090             zzlev(:,nlayer+1)=zzlay(:,nlayer)+(zzlay(:,nlayer)-zzlev(:,nlayer)) ! JVO 19 : We assume zzlev isn't reused later on
     1098
     1099            zzlev(:,nlayer+1)=zzlay(:,nlayer)+(zzlay(:,nlayer)-zzlev(:,nlayer)) ! JVO 19 : We assume zzlev isn't reused later on (could be done cleaner)
     1100
    10911101#ifdef USE_QTEST
    1092             dtpq(:,:,:) = 0.0 ! we want tpq to go only through mufi
     1102            dtpq(:,:,:) = 0.D0 ! we want tpq to go only through mufi
    10931103            call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,tpq,dtpq,zdqmufi)
    10941104            tpq(:,:,:) = tpq(:,:,:) + zdqmufi(:,:,:)*ptimestep ! only manipulation of tpq->*ptimestep here
    10951105#else
     1106           
    10961107            call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,zdqmufi)
     1108
    10971109            pdq(:,:,:) = pdq(:,:,:) + zdqmufi(:,:,:)
     1110
     1111            ! Sanity check ( way safer to be done here rather than within YAMMS )
     1112            ! Important : the sanity check intentionally include the former processes tendency !
     1113            ! NB : Despite this sanity check there might be still some unphysical values going through :
     1114            !        - Negatives, but harmless as it will be only for the output files
     1115            !          just remove them in post-proc.
     1116            !        - Weird unphysical ratio of m0 and m3, ok for now, but take care of them if
     1117            !          you want to compute optics from radii.
     1118            WHERE ( (pq(:,:,1)+pdq(:,:,1)*ptimestep < 0.D0) .OR. (pq(:,:,2)+pdq(:,:,2)*ptimestep < 0.D0) )
     1119                    pdq(:,:,1) = (epsilon(1.0)-1.D0)*pq(:,:,1)/ptimestep
     1120                    pdq(:,:,2) = (epsilon(1.0)-1.D0)*pq(:,:,2)/ptimestep
     1121            ENDWHERE
     1122            WHERE ( (pq(:,:,3)+pdq(:,:,3)*ptimestep < 0.D0) .OR. (pq(:,:,4)+pdq(:,:,4)*ptimestep < 0.D0) )
     1123                    pdq(:,:,3) = (epsilon(1.0)-1.D0)*pq(:,:,3)/ptimestep
     1124                    pdq(:,:,4) = (epsilon(1.0)-1.D0)*pq(:,:,4)/ptimestep
     1125            ENDWHERE
    10981126#endif
    10991127
    … …  
    11031131              call calmufi(ptimestep,zplevbar,zzlevbar,zplaybar,zzlaybar, &
    11041132                           gzlat,ztfibar,zqfibar,zdqfibar,zdqmufibar)
     1133              ! TODO : Add a sanity check here !
    11051134            endif
    11061135
    … …  
    12151244            ! v. Updates and positivity check
    12161245            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1217             zdqchi(:,:,:)  = 0.D0 ! -> for the microphysical tracers we need to update to 0 each time ;-)
     1246            zdqchi(:,:,:)  = 0.D0 ! -> dycchi is saved but for the nmicro tracers we must update to 0 at each step
    12181247
    12191248            do iq=1,nkim
    … …  
    13151344
    13161345      if(firstcall)then
    1317          zdtdyn(:,:)=0.0
    1318          zdudyn(:,:)=0.0
     1346         zdtdyn(:,:)=0.D0
     1347         zdudyn(:,:)=0.D0
    13191348      endif
    13201349
    … …  
    16741703        CALL send_xios_field("dqcond_C3H6",dyccond(:,:,17+nmicro))
    16751704        CALL send_xios_field("dqcond_C4H4",dyccond(:,:,21+nmicro))
    1676         CALL send_xios_field("dqcond_CH3CCH",dyccond(:,:,24+nmicro))
    1677         CALL send_xios_field("dqcond_C3H8",dyccond(:,:,25+nmicro))
    1678         CALL send_xios_field("dqcond_C4H2",dyccond(:,:,26+nmicro))
    1679         CALL send_xios_field("dqcond_C4H6",dyccond(:,:,27+nmicro))
    1680         CALL send_xios_field("dqcond_C4H10",dyccond(:,:,28+nmicro))
    1681         CALL send_xios_field("dqcond_AC6H6",dyccond(:,:,29+nmicro))
    1682         CALL send_xios_field("dqcond_HCN",dyccond(:,:,36+nmicro))
    1683         CALL send_xios_field("dqcond_CH3CN",dyccond(:,:,40+nmicro))
    1684         CALL send_xios_field("dqcond_HC3N",dyccond(:,:,42+nmicro))
    1685         CALL send_xios_field("dqcond_NCCN",dyccond(:,:,43+nmicro))
    1686         CALL send_xios_field("dqcond_C4N2",dyccond(:,:,44+nmicro))
     1705        CALL send_xios_field("dqcond_CH3CCH",dyccond(:,:,23+nmicro))
     1706        CALL send_xios_field("dqcond_C3H8",dyccond(:,:,24+nmicro))
     1707        CALL send_xios_field("dqcond_C4H2",dyccond(:,:,25+nmicro))
     1708        CALL send_xios_field("dqcond_C4H6",dyccond(:,:,26+nmicro))
     1709        CALL send_xios_field("dqcond_C4H10",dyccond(:,:,27+nmicro))
     1710        CALL send_xios_field("dqcond_AC6H6",dyccond(:,:,28+nmicro))
     1711        CALL send_xios_field("dqcond_HCN",dyccond(:,:,35+nmicro))
     1712        CALL send_xios_field("dqcond_CH3CN",dyccond(:,:,39+nmicro))
     1713        CALL send_xios_field("dqcond_HC3N",dyccond(:,:,41+nmicro))
     1714        CALL send_xios_field("dqcond_NCCN",dyccond(:,:,42+nmicro))
     1715        CALL send_xios_field("dqcond_C4N2",dyccond(:,:,43+nmicro))
    16871716
    16881717        ! Pseudo-evaporation flux (mol/mol/s)
Note: See TracChangeset for help on using the changeset viewer.