Changeset 4982


Ignore:
Timestamp:
Jun 15, 2024, 5:17:08 PM (4 days ago)
Author:
crisi
Message:

suppress isotope_params.def + update physiq_mod + proof of concept of 3rd dimension with reevap routine

Location:
LMDZ6/trunk/libf/phylmdiso
Files:
2 added
9 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmdiso/add_wake_tend.F90

    r4143 r4982  
    1 SUBROUTINE add_wake_tend(zddeltat, zddeltaq, zds, zddensaw, zddensw, zoccur, text, abortphy &
     1SUBROUTINE add_wake_tend(zddeltat, zddeltaq, zds, zdas, zddensw, zddensaw, zoccur, text, abortphy &
    22#ifdef ISO
    33        , zddeltaxt &
    … …  
    1313
    1414USE dimphy, ONLY: klon, klev
    15 USE phys_state_var_mod, ONLY: wake_deltat, wake_deltaq, wake_s, &
    16                               awake_dens, wake_dens
     15USE phys_state_var_mod, ONLY: wake_deltat, wake_deltaq, wake_s, awake_s, &
     16                              wake_dens, awake_dens
    1717
    1818USE print_control_mod, ONLY: prt_level
    … …  
    2626!------------
    2727  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: zddeltat, zddeltaq
    28   REAL, DIMENSION(klon),         INTENT (IN)         :: zds, zddensaw, zddensw
     28  REAL, DIMENSION(klon),         INTENT (IN)         :: zds, zdas, zddensw, zddensaw
    2929  INTEGER, DIMENSION(klon),      INTENT (IN)         :: zoccur
    3030  CHARACTER*(*),                 INTENT (IN)         :: text
    … …  
    7979           IF (zoccur(i) .GE. 1) THEN
    8080             wake_s(i)     = wake_s(i)    + zds(i)
     81             awake_s(i)    = awake_s(i)    + zdas(i)
     82             wake_dens(i)  = wake_dens(i) + zddensw(i)
    8183             awake_dens(i) = awake_dens(i) + zddensaw(i)
    82              wake_dens(i)  = wake_dens(i) + zddensw(i)
    8384           ELSE
    8485             wake_s(i)     = 0.
     86             awake_s(i)    = 0.
     87             wake_dens(i)  = 0.
    8588             awake_dens(i) = 0.
    86              wake_dens(i)  = 0.
    8789           ENDIF   ! (zoccur(i) .GE. 1)
    8890         END DO

  • LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90

    r4491 r4982  
    2020
    2121   !--- Variables not depending on isotopes
    22    REAL,    SAVE :: pxtmelt, pxtice, pxtmin, pxtmax
    23 !$OMP THREADPRIVATE(pxtmelt, pxtice, pxtmin, pxtmax)
    24    REAL,    SAVE :: tdifexp, tv0cin, thumxt1
    25 !$OMP THREADPRIVATE(tdifexp, tv0cin, thumxt1)
     22   REAL,    SAVE :: thumxt1
     23   !$OMP THREADPRIVATE(thumxt1)
    2624   INTEGER, SAVE :: ntot
    2725!$OMP THREADPRIVATE(ntot)
    … …  
    3028   REAL,    SAVE :: P_veg
    3129!$OMP THREADPRIVATE(P_veg)
    32    REAL,    SAVE :: musi, lambda_sursat
    33 !$OMP THREADPRIVATE(musi, lambda_sursat)
    34    REAL,    SAVE :: Kd
    35 !$OMP THREADPRIVATE(Kd)
    36    REAL,    SAVE :: rh_cste_surf_cond, T_cste_surf_cond
    37 !$OMP THREADPRIVATE(rh_cste_surf_cond, T_cste_surf_cond)
     30   REAL,    SAVE :: lambda_sursat
     31!$OMP THREADPRIVATE(lambda_sursat)
    3832   LOGICAL, SAVE :: bidouille_anti_divergence    ! T: regularly, xteau <- q to avoid slow drifts
    3933!$OMP THREADPRIVATE(bidouille_anti_divergence)
    … …  
    5448   REAL,    SAVE :: sstlatcrit, dsstlatcrit
    5549!$OMP THREADPRIVATE(sstlatcrit, dsstlatcrit)
    56    REAL,    SAVE :: deltaO18_oce
    57 !$OMP THREADPRIVATE(deltaO18_oce)
    5850   INTEGER, SAVE :: albedo_prescrit              ! 0: default ; 1: constant albedo
    5951!$OMP THREADPRIVATE(albedo_prescrit)
    … …  
    8880   REAL,    SAVE :: fac_modif_evaoce
    8981!$OMP THREADPRIVATE(fac_modif_evaoce)
     82   REAL,    SAVE :: deltaO18_oce
     83!$OMP THREADPRIVATE(deltaO18_oce)
    9084   INTEGER, SAVE :: ok_bidouille_wake
    9185!$OMP THREADPRIVATE(ok_bidouille_wake)
    … …  
    106100                    alpha_liq_sol, Rdefault, Rmethox
    107101!$OMP THREADPRIVATE(alpha_liq_sol, Rdefault, Rmethox)
    108    REAL, SAVE ::    fac_coeff_eq17_liq, fac_coeff_eq17_ice
    109 !$OMP THREADPRIVATE(fac_coeff_eq17_liq, fac_coeff_eq17_ice)
     102!   REAL, SAVE ::    fac_coeff_eq17_liq, fac_coeff_eq17_ice
     103!!$OMP THREADPRIVATE(fac_coeff_eq17_liq, fac_coeff_eq17_ice)
     104
     105   !--- H2[18]O reference
     106   REAL, PARAMETER :: fac_enrichoce18=0.0005
     107   REAL, PARAMETER :: alpha_liq_sol_O18=1.00291
     108   REAL, PARAMETER :: talph1_O18=1137.
     109   REAL, PARAMETER :: talph2_O18=-0.4156
     110   REAL, PARAMETER :: talph3_O18=-2.0667E-3
     111   REAL, PARAMETER :: talps1_O18=11.839
     112   REAL, PARAMETER :: talps2_O18=-0.028244
     113   REAL, PARAMETER :: tdifrel_O18=1./0.9723
     114   REAL, PARAMETER :: tkcin0_O18=0.006
     115   REAL, PARAMETER :: tkcin1_O18=0.000285
     116   REAL, PARAMETER :: tkcin2_O18=0.00082
     117   REAL, PARAMETER :: fac_coeff_eq17_liq=0.529
     118   REAL, PARAMETER :: fac_coeff_eq17_ice=0.529
     119
     120   !---- Parameters that do not depend on the nature of water isotopes:
     121   REAL, PARAMETER :: pxtmelt = 273.15 ! temperature at which ice formation starts
     122   REAL, PARAMETER :: pxtice  = 273.15-10.0 ! -- temperature at which all condensate is ice:
     123   REAL, PARAMETER :: pxtmin = 273.15 - 120.0   ! On ne calcule qu'au dessus de -120°C
     124   REAL, PARAMETER :: pxtmax = 273.15 +  60.0   ! On ne calcule qu'au dessus de +60°C
     125   REAL, PARAMETER :: tdifexp = 0.58 ! -- a constant for alpha_eff for equilibrium below cloud base:
     126   REAL, PARAMETER :: tv0cin  = 7.0 ! wind threshold (m/s) for smooth/rough regime (Merlivat and Jouzel 1979)
     127   REAL, PARAMETER :: musi=1.0  ! facteurs lambda et mu dans Si=musi-lambda*T
     128   REAL, PARAMETER :: Kd=2.5e-9 ! m2/s ! diffusion dans le sol
     129   REAL, PARAMETER :: rh_cste_surf_cond = 0.6 ! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir
     130   REAL, PARAMETER :: T_cste_surf_cond = 288.0
     131
    110132
    111133   !--- Negligible lower thresholds: no need to check for absurd values under these lower limits
    … …  
    140162   INTEGER :: ixt
    141163
    142    !--- H2[18]O reference
    143    REAL :: fac_enrichoce18, alpha_liq_sol_O18, &
    144            talph1_O18, talph2_O18, talph3_O18, talps1_O18, talps2_O18, &
    145            tkcin0_O18, tkcin1_O18, tkcin2_O18, tdifrel_O18 
    146 
     164 
    147165   !--- For H2[17]O
    148166   REAL    :: fac_kcin, pente_MWL
    … …  
    152170   LOGICAL, PARAMETER ::   ok_nocinsat = .FALSE. ! if T: no sursaturation effect for ice
    153171   LOGICAL, PARAMETER :: Rdefault_smow = .FALSE. ! if T: Rdefault=smow; if F: nul
     172   LOGICAL, PARAMETER :: tnat1 = .TRUE. ! If T: all tnats are 1.
    154173
    155174   !--- For [3]H
    … …  
    157176
    158177   CHARACTER(LEN=maxlen) :: modname, sxt
    159    REAL, ALLOCATABLE :: tmp(:)
    160178
    161179   modname = 'iso_init'
    … …  
    214232      CALL get_in('lat_max_albedo', lat_max_albedo,  100.)
    215233   END IF
    216    deltaO18_oce=0.0
    217234   CALL get_in('deltaP_BL',           deltaP_BL,     10.0)
    218235   CALL get_in('ruissellement_pluie', ruissellement_pluie, 0)
    … …  
    249266   CALL get_in('ok_prod_nucl_tritium', ok_prod_nucl_tritium, .FALSE., .FALSE.)
    250267
    251    !--------------------------------------------------------------
    252    ! Parameters that do not depend on the nature of water isotopes:
    253    !--------------------------------------------------------------
    254    ! -- temperature at which ice condensate starts to form (valeur ECHAM?):
    255    pxtmelt = 273.15
    256 
    257    ! -- temperature at which all condensate is ice:
    258    pxtice  = 273.15-10.0
    259 
    260    !- -- test PHASE
    261 !   pxtmelt = 273.15 - 10.0
    262 !   pxtice  = 273.15 - 30.0
    263 
    264    ! -- minimum temperature to calculate fractionation coeff
    265    pxtmin = 273.15 - 120.0   ! On ne calcule qu'au dessus de -120°C
    266    pxtmax = 273.15 +  60.0   ! On ne calcule qu'au dessus de +60°C
    267    !    Remarque: les coeffs ont ete mesures seulement jusq'à -40!
    268 
    269    ! -- a constant for alpha_eff for equilibrium below cloud base:
    270    tdifexp = 0.58
    271    tv0cin  = 7.0
    272 
    273    ! facteurs lambda et mu dans Si=musi-lambda*T
    274    musi=1.0
    275    if (ok_nocinsat) lambda_sursat = 0.0          ! no sursaturation effect
    276 
    277    ! diffusion dans le sol
    278    Kd=2.5e-9 ! m2/s   
    279 
    280    ! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir
    281    rh_cste_surf_cond = 0.6
    282     T_cste_surf_cond = 288.0
     268   ! Ocean composition
     269   CALL get_in('deltaO18_oce',  deltaO18_oce, 0.0)
    283270   
    284271   CALL msg('iso_O18, iso_HDO, iso_eau = '//TRIM(strStack(int2str([iso_O18, iso_HDO, iso_eau]))), modname)
    285272
    286273   !--------------------------------------------------------------
    287    ! Parameters that depend on the nature of water isotopes:
     274   ! Isotope fractionation factors and a few isotopic constants
    288275   !--------------------------------------------------------------
    289    IF(getKey('tnat',    tnat,    isoName)) CALL abort_physic(modname, 'can''t get tnat',    1)
    290    IF(getKey('toce',    toce,    isoName)) CALL abort_physic(modname, 'can''t get toce',    1)
    291    IF(getKey('tcorr',   tcorr,   isoName)) CALL abort_physic(modname, 'can''t get tcorr',   1)
    292    IF(getKey('talph1',  talph1,  isoName)) CALL abort_physic(modname, 'can''t get talph1',  1)
    293    IF(getKey('talph2',  talph2,  isoName)) CALL abort_physic(modname, 'can''t get talph2',  1)
    294    IF(getKey('talph3',  talph3,  isoName)) CALL abort_physic(modname, 'can''t get talph3',  1)
    295    IF(getKey('talps1',  talps1,  isoName)) CALL abort_physic(modname, 'can''t get talps1',  1)
    296    IF(getKey('talps2',  talps2,  isoName)) CALL abort_physic(modname, 'can''t get talps2',  1)
    297    IF(getKey('tkcin0',  tkcin0,  isoName)) CALL abort_physic(modname, 'can''t get tkcin0',  1)
    298    IF(getKey('tkcin1',  tkcin1,  isoName)) CALL abort_physic(modname, 'can''t get tkcin1',  1)
    299    IF(getKey('tkcin2',  tkcin2,  isoName)) CALL abort_physic(modname, 'can''t get tkcin2',  1)
    300    IF(getKey('tdifrel', tdifrel, isoName)) CALL abort_physic(modname, 'can''t get tdifrel', 1)
    301    IF(getKey('alpha_liq_sol', alpha_liq_sol, isoName)) CALL abort_physic(modname, 'can''t get alpha_liq_sol',  1)
    302    IF(getKey('Rdefault',Rdefault,isoName)) CALL abort_physic(modname, 'can''t get Rdefault',1)
    303    IF(getKey('Rmethox', Rmethox, isoName)) CALL abort_physic(modname, 'can''t get Rmethox', 1)
     276   ALLOCATE(tkcin0(niso))
     277   ALLOCATE(tkcin1(niso))
     278   ALLOCATE(tkcin2(niso))
     279   ALLOCATE(tnat(niso))
     280   ALLOCATE(tdifrel(niso))
     281   ALLOCATE(toce(niso))
     282   ALLOCATE(tcorr(niso))
     283   ALLOCATE(talph1(niso))
     284   ALLOCATE(talph2(niso))
     285   ALLOCATE(talph3(niso))
     286   ALLOCATE(talps1(niso))
     287   ALLOCATE(talps2(niso))
     288   ALLOCATE(alpha_liq_sol(niso))
     289   ALLOCATE(Rdefault(niso))
     290   ALLOCATE(Rmethox(niso))
     291
     292   do ixt=1,niso
     293     if (ixt.eq.iso_HTO) then  ! Tritium
     294       tkcin0(ixt) = 0.01056
     295       tkcin1(ixt) = 0.0005016
     296       tkcin2(ixt) = 0.0014432
     297       if (tnat1) then
     298               tnat(ixt)=1
     299       else
     300               tnat(ixt)=0.
     301       endif
     302       toce(ixt)=4.0E-19 ! rapport T/H = 0.2 TU Dreisigacker and Roether 1978
     303       tcorr(ixt)=1.
     304       tdifrel(ixt)=1./0.968
     305       talph1(ixt)=46480.
     306       talph2(ixt)=-103.87
     307       talph3(ixt)=0.
     308       talps1(ixt)=46480.
     309       talps2(ixt)=-103.87
     310       alpha_liq_sol(ixt)=1.
     311       Rmethox(ixt)=0.0
     312     endif
     313     if (ixt.eq.iso_O17) then  ! O17
     314       pente_MWL=0.528
     315       tdifrel(ixt)=1./0.98555 ! valeur utilisée en 1D et dans modèle de LdG ! tdifrel(ixt)=1./0.985452 ! donné par Amaelle
     316       fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0) ! fac_kcin=0.5145 ! donné par Amaelle
     317       tkcin0(ixt) = tkcin0_O18*fac_kcin
     318       tkcin1(ixt) = tkcin1_O18*fac_kcin
     319       tkcin2(ixt) = tkcin2_O18*fac_kcin
     320       if (tnat1) then
     321               tnat(ixt)=1
     322       else
     323               tnat(ixt)=0.004/100. ! O17 représente 0.004% de l'oxygène
     324       endif
     325       toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)**pente_MWL
     326       tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL ! donné par Amaelle           
     327       talph1(ixt)=talph1_O18
     328       talph2(ixt)=talph2_O18
     329       talph3(ixt)=talph3_O18
     330       talps1(ixt)=talps1_O18
     331       talps2(ixt)=talps2_O18     
     332       alpha_liq_sol(ixt)=(alpha_liq_sol_O18)**fac_coeff_eq17_liq
     333       Rdefault(ixt)=tnat(ixt)*(-3.15/1000.0+1.0)
     334       Rmethox(ixt)=(230./1000.+1.)*tnat(ixt) !Zahn et al 2006
     335     endif
     336     if (ixt.eq.iso_O18) then  ! Oxygene18
     337       tkcin0(ixt) = tkcin0_O18
     338       tkcin1(ixt) = tkcin1_O18
     339       tkcin2(ixt) = tkcin2_O18
     340       if (tnat1) then
     341               tnat(ixt)=1
     342       else
     343               tnat(ixt)=2005.2E-6
     344       endif
     345       toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)
     346       tcorr(ixt)=1.0+fac_enrichoce18
     347       tdifrel(ixt)=tdifrel_O18
     348       talph1(ixt)=talph1_O18
     349       talph2(ixt)=talph2_O18
     350       talph3(ixt)=talph3_O18
     351       talps1(ixt)=talps1_O18
     352       talps2(ixt)=talps2_O18
     353       alpha_liq_sol(ixt)=alpha_liq_sol_O18   
     354       Rdefault(ixt)=tnat(ixt)*(-6.0/1000.0+1.0)
     355       Rmethox(ixt)=(130./1000.+1.)*tnat(ixt) !Zahn et al 2006 
     356     endif
     357     if (ixt.eq.iso_HDO) then ! Deuterium
     358       pente_MWL=8.0
     359       tdifrel(ixt)=1./0.9755 !          fac_kcin=0.88
     360       fac_kcin= (tdifrel(ixt)-1)/(tdifrel_O18-1)
     361       tkcin0(ixt) = tkcin0_O18*fac_kcin
     362       tkcin1(ixt) = tkcin1_O18*fac_kcin
     363       tkcin2(ixt) = tkcin2_O18*fac_kcin
     364       if (tnat1) then
     365               tnat(ixt)=1
     366       else
     367               tnat(ixt)=155.76E-6
     368       endif
     369       toce(ixt)=tnat(ixt)*(1.0+pente_MWL*deltaO18_oce/1000.0)
     370       tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL         
     371       talph1(ixt)=24844.
     372       talph2(ixt)=-76.248
     373       talph3(ixt)=52.612E-3
     374       talps1(ixt)=16288.
     375       talps2(ixt)=-0.0934
     376       !ZXalpha_liq_sol=1.0192 ! Weston, Ralph, 1955
     377       alpha_liq_sol(ixt)=1.0212
     378       ! valeur de Lehmann & Siegenthaler, 1991, Journal of
     379       ! Glaciology, vol 37, p 23
     380       Rdefault(ixt)=tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0)
     381       Rmethox(ixt)=tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006
     382     endif
     383     if (ixt.eq.iso_eau) then ! Oxygene16
     384       tkcin0(ixt) = 0.0
     385       tkcin1(ixt) = 0.0
     386       tkcin2(ixt) = 0.0
     387       tnat(ixt)=1.
     388       toce(ixt)=tnat(ixt)
     389       tcorr(ixt)=1.0
     390       tdifrel(ixt)=1.
     391       talph1(ixt)=0.
     392       talph2(ixt)=0.
     393       talph3(ixt)=0.
     394       talps1(ixt)=0.
     395       talph3(ixt)=0.
     396       alpha_liq_sol(ixt)=1.
     397       Rdefault(ixt)=tnat(ixt)*1.0
     398       Rmethox(ixt)=1.0
     399     endif
     400   enddo ! ixt=1,niso
    304401
    305402   IF(.NOT.Rdefault_smow) then
    … …  
    308405   ENDIF
    309406   write(*,*) 'Rdefault=',Rdefault
     407   write(*,*) 'toce=',toce
    310408
    311409   !--- Sensitivity test: no kinetic effect in sfc evaporation

  • LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90

    r4491 r4982  
    15251525#endif       
    15261526       pxtfra=max(min(pxtfra,alpha_max),0.0)
    1527 
    15281527
    15291528      end subroutine fractcalk_liq
    … …  
    1592215921
    1592315922      ! verif
    15924 !      text="phyisoetat0 67"
    15925 !      write(*,*) 'snow(8,1)=',snow(8,1)
    15926 !      write(*,*) 'xtsnow(4,8,1)=',xtsnow(4,8,1)
    1592715923#ifdef ISOVERIF
    1592815924      do i=1,klon
    … …  
    1593415930         enddo !do ixt=1,niso
    1593515931      enddo !do i=1,klon
    15936 #endif     
    15937 #ifdef ISOVERIF
    1593815932      do i=1,klon
    1593915933         if (iso_eau.gt.0) then
    … …  
    1602116015         endif
    1602216016       enddo !do i=1,klon
    16023 
    1602416017#endif
    1602516018      !end verif
    … …  
    1612816121          deltaD_run_off_lic_0(ixt)=deltaD_sol(ixt)
    1612916122          deltaD_land_ice(ixt)=deltaD_snow(ixt)
    16130           call fractcalk_liq(ixt, 283.0, alpha(ixt))           
     16123          call fractcalk_liq(ixt, 283.0, alpha(ixt))   
    1613116124        enddo !do ixt=1,niso
    1613216125        call calcul_kcin(2.0,kcin)
    … …  
    1883018823        if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO)) then
    1883118824            if (q.gt.ridicule) then
     18825                    write(*,*) 'xt,q=',xt,q
     18826                    write(*,*) 'alpha=',alpha
     18827                    write(*,*) 'toce,kcin,h0=',toce,kcin,h0
     18828                    write(*,*) 'RMerlivat=',RMerlivat
    1883218829                call iso_verif_aberrant_encadre( xt/q, 'isotopes_routines_mod 18930b: iso_init_ideal')
    1883318830            endif
    … …  
    1890218899end subroutine appel_stewart_debug
    1890318900
     18901
     18902subroutine dispatch(klon,klev,qx,q_seri,xt_seri,ql_seri,xtl_seri,qs_seri,xts_seri)
     18903
     18904use infotrac_phy, ONLY: nqtot,nqo,ivap,iliq,isol,iqIsoPha,ntraciso=>ntiso
     18905implicit none
     18906
     18907! inputs
     18908integer, intent(in) :: klon,klev
     18909real,dimension(klon,klev,nqtot), intent(in) ::qx
     18910
     18911! outputs
     18912real,dimension(klon,klev), intent(out) ::q_seri,ql_seri,qs_seri
     18913real,dimension(ntraciso,klon,klev), intent(out) :: xt_seri,xtl_seri,xts_seri
     18914
     18915! locals
     18916integer :: i,k,ixt
     18917
     18918do k=1,klev
     18919do i=1,klon
     18920    q_seri(i,k)  = qx(i,k,ivap)
     18921    ql_seri(i,k) = qx(i,k,iliq)
     18922    IF (nqo.EQ.2) THEN             !--vapour and liquid only
     18923             qs_seri(i,k) = 0.
     18924    ELSE IF (nqo.ge.3) THEN        !--vapour, liquid and ice
     18925             qs_seri(i,k) = qx(i,k,isol)
     18926    ENDIF
     18927    do ixt=1,ntraciso
     18928          xt_seri(ixt,i,k)  = qx(i,k,iqIsoPha(ixt,ivap))
     18929          xtl_seri(ixt,i,k) = qx(i,k,iqIsoPha(ixt,iliq))
     18930          if (nqo.eq.2) then
     18931             xts_seri(ixt,i,k) = 0.
     18932          else if (nqo.eq.3) then
     18933             xts_seri(ixt,i,k) = qx(i,k,iqIsoPha(ixt,isol))
     18934          endif
     18935    enddo !do ixt=1,niso
     18936
     18937enddo
     18938enddo
     18939
     18940end subroutine dispatch
     18941
     18942subroutine together(klon,klev,qx,q_seri,xt_seri,ql_seri,xtl_seri,qs_seri,xts_seri)
     18943
     18944use infotrac_phy, ONLY: nqtot,nqo,ivap,iliq,isol,iqIsoPha,ntraciso=>ntiso
     18945implicit none
     18946
     18947! inputs
     18948integer, intent(in) :: klon,klev
     18949real,dimension(klon,klev), intent(in) ::q_seri,ql_seri,qs_seri
     18950real,dimension(ntraciso,klon,klev), intent(in) :: xt_seri,xtl_seri,xts_seri
     18951
     18952! inputs
     18953real,dimension(klon,klev,nqtot), intent(out) ::qx
     18954
     18955! locals
     18956integer :: i,k,ixt
     18957
     18958do k=1,klev
     18959do i=1,klon 
     18960    qx(i,k,ivap)  = q_seri(i,k)
     18961    qx(i,k,iliq) = ql_seri(i,k)
     18962    IF (nqo.ge.3) THEN        !--vapour, liquid and ice
     18963             qx(i,k,isol) = qs_seri(i,k)
     18964    ENDIF
     18965    do ixt=1,ntraciso
     18966          qx(i,k,iqIsoPha(ixt,ivap)) = xt_seri(ixt,i,k) 
     18967          qx(i,k,iqIsoPha(ixt,iliq)) = xtl_seri(ixt,i,k)
     18968          if (nqo.ge.3) then
     18969             qx(i,k,iqIsoPha(ixt,isol)) = xts_seri(ixt,i,k)
     18970          endif
     18971    enddo !do ixt=1,niso
     18972
     18973enddo
     18974enddo
     18975
     18976end subroutine together
     18977
     18978
    1890418979END MODULE isotopes_routines_mod
    1890518980#endif

  • LMDZ6/trunk/libf/phylmdiso/isotopes_verif_mod.F90

    r4491 r4982  
    10421042            write(*,*) 'deltaD=',deltaD
    10431043            write(*,*) 'Dexcess=',dexcess
     1044            write(*,*) 'tnat=',tnat
    10441045!            stop
    10451046            iso_verif_O18_aberrant_nostop=1

  • LMDZ6/trunk/libf/phylmdiso/isotrac_routines_mod.F90

    r4491 r4982  
    681681              Eqi_prime_cas(il)=Eqi_prime(cas(il)) &
    682682     &           *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
    683               Eqi_cas(il)=Eqi(il) &
     683              Eqi_cas(il)=Eqi(cas(il)) & ! corr bug Camille 15 juin 2024
    684684     &           *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
    685685            else

  • LMDZ6/trunk/libf/phylmdiso/phyetat0_mod.F90

    r4619 r4982  
    2929       solsw, solswfdiff, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
    3030       wake_deltat, wake_delta_pbl_TKE, delta_tsurf, beta_aridity, wake_fip, wake_pe, &
    31        wake_s, wake_dens, awake_dens, cv_gen, zgam, zmax0, zmea, zpic, zsig, &
     31       wake_s, awake_s, wake_dens, awake_dens, cv_gen, zgam, zmax0, zmea, zpic, zsig, &
    3232#ifdef ISO
    3333       fxtevap, xtsol, xt_ancien, xtl_ancien, xts_ancien, wake_deltaxt, &
    … …  
    4848  USE time_phylmdz_mod, ONLY: init_iteration, pdtphys, itau_phy
    4949  USE wxios, ONLY: missing_val_xios => missing_val, using_xios
    50   use netcdf, only: nf90_fill_real
     50  use netcdf, only: missing_val_netcdf => nf90_fill_real
    5151  use config_ocean_skin_m, only: activate_ocean_skin
    5252#ifdef ISO
    … …  
    111111  REAL Rland_ice(niso,klon)
    112112#endif
     113
     114  IF (using_xios) THEN
     115    missing_val=missing_val_xios
     116  ELSE
     117    missing_val=missing_val_netcdf
     118  ENDIF
     119
    113120  ! FH1D
    114121  !     real iolat(jjm+1)
    … …  
    116123
    117124  ! Ouvrir le fichier contenant l'etat initial:
    118   IF (using_xios) THEN
    119     missing_val = missing_val_xios
    120   ELSE
    121     missing_val =  nf90_fill_real
    122   ENDIF
    123125
    124126  CALL open_startphy(fichnom)
    … …  
    323325
    324326!===================================================================
     327! Lecture dans le cas iflag_pbl_surface =1
     328!===================================================================
     329
     330   if ( iflag_physiq <= 1 ) then
     331!===================================================================
    325332  ! Lecture des temperatures du sol profond:
    326333!===================================================================
    … …  
    350357  found=phyetat0_get(snow_fall,"snow_f","snow fall",0.)
    351358  found=phyetat0_get(rain_fall,"rain_f","rain fall",0.)
    352 
    353359  IF (ok_bs) THEN
    354360     found=phyetat0_get(bs_fall,"bs_f","blowing snow fall",0.)
    … …  
    404410  ENDIF
    405411
     412  endif ! iflag_physiq <= 1
     413
    406414  ! Lecture de l'age de la neige:
    407415  found=phyetat0_srf(agesno,"AGESNO","SNOW AGE",0.001)
    … …  
    427435     prbsw_ancien(:)=0.
    428436  ENDIF
    429 
    430437
    431438  ! Ehouarn: addtional tests to check if t_ancien, q_ancien contain
    … …  
    450457  ENDIF
    451458
    452 
    453459  found=phyetat0_get(clwcon,"CLWCON","CLWCON",0.)
    454460  found=phyetat0_get(rnebcon,"RNEBCON","RNEBCON",0.)
    … …  
    485491  found=phyetat0_get(wake_deltaq,"WAKE_DELTAQ","Delta hum. wake/env",0.)
    486492  found=phyetat0_get(wake_s,"WAKE_S","Wake frac. area",0.)
     493  found=phyetat0_get(awake_s,"AWAKE_S","Active Wake frac. area",0.)
    487494!jyg<
    488495!  Set wake_dens to -1000. when there is no restart so that the actual
    … …  
    664671!        write(*,*) 'xtsnow(:,994,2)=',xtsnow(:,994,2)
    665672!#endif
    666 
     673  if ( iflag_physiq <= 1 ) then
    667674  CALL pbl_surface_init(fder, snow, qsurf, tsoil)
    668675#ifdef ISO
    669676  CALL pbl_surface_init_iso(xtsnow,Rland_ice)
    670677#endif
     678  endif
    671679
    672680  ! Initialize module ocean_cpl_mod for the case of coupled ocean

  • LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90

    r4920 r4982  
    1414      REAL, SAVE, ALLOCATABLE :: ql_seri(:,:),qs_seri(:,:)
    1515      !$OMP THREADPRIVATE(ql_seri,qs_seri)
     16      REAL, SAVE, ALLOCATABLE :: qx_seri(:,:,:)
     17      !$OMP THREADPRIVATE(qx_seri)
    1618      REAL, SAVE, ALLOCATABLE :: qbs_seri(:,:)
    1719      !$OMP THREADPRIVATE(qbs_seri)
    … …  
    6466      REAL, SAVE, ALLOCATABLE :: d_t_eva(:,:),d_q_eva(:,:),d_ql_eva(:,:),d_qi_eva(:,:)
    6567      !$OMP THREADPRIVATE(d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
     68      REAL, SAVE, ALLOCATABLE :: d_qx_eva(:,:,:)
     69      !$OMP THREADPRIVATE(d_qx_eva)
    6670      REAL, SAVE, ALLOCATABLE :: d_t_lscst(:,:),d_q_lscst(:,:)
    6771      !$OMP THREADPRIVATE(d_t_lscst,d_q_lscst)
    … …  
    124128      REAL, SAVE, ALLOCATABLE :: xts_seri(:,:,:)
    125129      !$OMP THREADPRIVATE( xts_seri)
     130      REAL, SAVE, ALLOCATABLE :: xtbs_seri(:,:,:)
     131      !$OMP THREADPRIVATE( xtbs_seri)
    126132      REAL, SAVE, ALLOCATABLE :: d_xt_eva(:,:,:)
    127133      !$OMP THREADPRIVATE( d_xt_eva)
    … …  
    134140      REAL, SAVE, ALLOCATABLE :: d_xt_dyn(:,:,:)
    135141      !$OMP THREADPRIVATE( d_xt_dyn)
    136       REAL, SAVE, ALLOCATABLE :: d_xtl_dyn(:,:,:), d_xts_dyn(:,:,:)
    137       !$OMP THREADPRIVATE(d_xtl_dyn, d_xts_dyn)
     142      REAL, SAVE, ALLOCATABLE :: d_xtl_dyn(:,:,:), d_xts_dyn(:,:,:), d_xtbs_dyn(:,:,:)
     143      !$OMP THREADPRIVATE(d_xtl_dyn, d_xts_dyn, d_xtbs_dyn)
    138144      REAL, SAVE, ALLOCATABLE :: d_xt_con(:,:,:)
    139145      !$OMP THREADPRIVATE( d_xt_con)
    … …  
    292298!$OMP THREADPRIVATE(toplwad0_aerop, sollwad0_aerop)
    293299
     300!AI 08 2023 ajout pour Ecrad
     301      REAL,ALLOCATABLE,SAVE :: topswad_aero_s2(:), solswad_aero_s2(:)
     302!$OMP THREADPRIVATE(topswad_aero_s2, solswad_aero_s2)
     303      REAL,ALLOCATABLE,SAVE :: topswai_aero_s2(:), solswai_aero_s2(:)
     304!$OMP THREADPRIVATE(topswai_aero_s2, solswai_aero_s2)
     305      REAL,ALLOCATABLE,SAVE :: topswad0_aero_s2(:), solswad0_aero_s2(:)
     306!$OMP THREADPRIVATE(topswad0_aero_s2, solswad0_aero_s2)
     307      REAL,ALLOCATABLE,SAVE :: topsw_aero_s2(:,:), topsw0_aero_s2(:,:)
     308!$OMP THREADPRIVATE(topsw_aero_s2, topsw0_aero_s2)
     309      REAL,ALLOCATABLE,SAVE :: solsw_aero_s2(:,:), solsw0_aero_s2(:,:)
     310!$OMP THREADPRIVATE(solsw_aero_s2, solsw0_aero_s2)
     311      REAL,ALLOCATABLE,SAVE :: topswcf_aero_s2(:,:), solswcf_aero_s2(:,:)
     312!$OMP THREADPRIVATE(topswcf_aero_s2, solswcf_aero_s2)
     313! additional LW variables CK
     314      REAL,ALLOCATABLE,SAVE :: toplwad_aero_s2(:), sollwad_aero_s2(:)
     315!$OMP THREADPRIVATE(toplwad_aero_s2, sollwad_aero_s2)
     316      REAL,ALLOCATABLE,SAVE :: toplwai_aero_s2(:), sollwai_aero_s2(:)
     317!$OMP THREADPRIVATE(toplwai_aero_s2, sollwai_aero_s2)
     318      REAL,ALLOCATABLE,SAVE :: toplwad0_aero_s2(:), sollwad0_aero_s2(:)
     319!$OMP THREADPRIVATE(toplwad0_aero_s2, sollwad0_aero_s2)
     320
    294321!Ajout de celles n??cessaires au phys_output_write_mod
    295322      REAL, SAVE, ALLOCATABLE :: tal1(:), pal1(:), pab1(:), pab2(:)
    … …  
    300327!$OMP THREADPRIVATE(sens, flwp, fiwp)
    301328!!
    302 !FC
     329!FC 
    303330      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: zxfluxt, zxfluxq
    304331!$OMP THREADPRIVATE(zxfluxt, zxfluxq)
    … …  
    348375!$OMP THREADPRIVATE(cdragm, cdragh)
    349376      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: cldh, cldl, cldm, cldq, cldt, qsat2m
    350 !$OMP THREADPRIVATE(cldh, cldl, cldm, cldq, cldt, qsat2m )
     377!$OMP THREADPRIVATE(cldh, cldl, cldm, cldq, cldt, qsat2m)
    351378!AS: cldhjn, cldljn, cldmjn,cldtjn pas utilisés en tant que variables, juste noms de diagnostics
    352379      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: JrNt
    … …  
    573600!$OMP THREADPRIVATE(phiwriteSTD, qwriteSTD, twriteSTD, rhwriteSTD)
    574601
     602
    575603! ug et d'autres encore:
    576604      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: beta_prec
    … …  
    580608      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: pfraclr,pfracld
    581609!$OMP THREADPRIVATE(pfraclr,pfracld)
    582       REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: distcltop
    583 !$OMP THREADPRIVATE(distcltop)
    584       REAL, SAVE, ALLOCATABLE :: temp_cltop(:,:)
    585 !$OMP THREADPRIVATE(temp_cltop)
    586610
    587611! variables de sorties MM
    … …  
    645669      REAL, SAVE, ALLOCATABLE :: fcontrP(:,:)
    646670      !$OMP THREADPRIVATE(fcontrP)
     671      REAL, SAVE, ALLOCATABLE :: distcltop(:,:)
     672      !$OMP THREADPRIVATE(distcltop)
     673      REAL, SAVE, ALLOCATABLE :: temp_cltop(:,:)
     674      !$OMP THREADPRIVATE(temp_cltop)
     675
    647676
    648677!--POPRECIP variables
    … …  
    675704
    676705
    677 
     706     
    678707
    679708
    … …  
    681710!
    682711! variables for stratospheric aerosol
     712      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: d_q_emiss
     713!$OMP THREADPRIVATE(d_q_emiss)
    683714      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: R2SO4
    684715!$OMP THREADPRIVATE(R2SO4)
    … …  
    695726      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: SO2_lifetime
    696727!$OMP THREADPRIVATE(SO2_lifetime)
     728      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: H2SO4_lifetime
     729!$OMP THREADPRIVATE(H2SO4_lifetime)
     730      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: O3_clim
     731!$OMP THREADPRIVATE(O3_clim)
    697732      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: alpha_bin
    698733!$OMP THREADPRIVATE(alpha_bin)
    … …  
    759794SUBROUTINE phys_local_var_init
    760795USE dimphy
    761 USE infotrac_phy, ONLY : nbtr
     796USE infotrac_phy, ONLY : nbtr,nqtot
    762797#ifdef ISO
    763798USE infotrac_phy, ONLY : ntraciso=>ntiso,niso
    … …  
    770805IMPLICIT NONE
    771806      ALLOCATE(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev), qbs_seri(klon,klev))
     807      ALLOCATE(qx_seri(klon,klev,nqtot))
    772808      ALLOCATE(u_seri(klon,klev),v_seri(klon,klev))
    773809      ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
     810      ALLOCATE(pbl_eps(klon,klev+1,nbsrf+1))
     811      pbl_eps(:,:,:)=0.
    774812      l_mix(:,:,:)=0.;l_mixmin(:,:,:)=0.;wprime(:,:,:)=0. ! doit etre initialse car pas toujours remplis
    775       ALLOCATE(pbl_eps(klon,klev+1,nbsrf+1))     
     813      ALLOCATE(rhcl(klon,klev))
    776814      ALLOCATE(tr_seri(klon,klev,nbtr))
    777815      ALLOCATE(d_t_dyn(klon,klev),d_q_dyn(klon,klev))
    … …  
    795833      ALLOCATE(d_u_ajs(klon,klev),d_v_ajs(klon,klev))
    796834      ALLOCATE(d_t_eva(klon,klev),d_q_eva(klon,klev))
     835      ALLOCATE(d_qx_eva(klon,klev,nqtot))
    797836      ALLOCATE(d_ql_eva(klon,klev),d_qi_eva(klon,klev))
    798837      ALLOCATE(d_t_lscst(klon,klev),d_q_lscst(klon,klev))
    … …  
    808847      allocate(xtl_seri(ntraciso,klon,klev))
    809848      allocate(xts_seri(ntraciso,klon,klev))
     849      allocate(xtbs_seri(ntraciso,klon,klev))
    810850      allocate(d_xt_dyn(ntraciso,klon,klev))
    811851      allocate(d_xtl_dyn(ntraciso,klon,klev))
    812852      allocate(d_xts_dyn(ntraciso,klon,klev))
     853      allocate(d_xtbs_dyn(ntraciso,klon,klev))
    813854      allocate(d_xt_con(ntraciso,klon,klev))
    814855      allocate(d_xt_wake(ntraciso,klon,klev))
    … …  
    841882      ALLOCATE(d_u_lif(klon,klev),d_v_lif(klon,klev))
    842883      ALLOCATE(d_ts(klon,nbsrf), d_tr(klon,klev,nbtr))
     884
     885! aerosols
     886      ALLOCATE(m_allaer(klon,klev,naero_tot))
    843887! Special RRTM
    844888      ALLOCATE(ZLWFT0_i(klon,klev+1),ZSWFT0_i(klon,klev+1),ZFLDN0(klon,klev+1))
    … …  
    919963      ALLOCATE(toplwad0_aerop(klon), sollwad0_aerop(klon))
    920964
     965!AI Ajout Ecrad (3Deffect)
     966      ALLOCATE(topswad_aero_s2(klon), solswad_aero_s2(klon))
     967      ALLOCATE(topswai_aero_s2(klon), solswai_aero_s2(klon))
     968      ALLOCATE(topswad0_aero_s2(klon), solswad0_aero_s2(klon))
     969      ALLOCATE(topsw_aero_s2(klon,naero_grp), topsw0_aero_s2(klon,naero_grp))
     970      ALLOCATE(solsw_aero_s2(klon,naero_grp), solsw0_aero_s2(klon,naero_grp))
     971      ALLOCATE(topswcf_aero_s2(klon,naero_grp), solswcf_aero_s2(klon,naero_grp))
     972! additional LW variables CK
     973      ALLOCATE(toplwad_aero_s2(klon), sollwad_aero_s2(klon))
     974      ALLOCATE(toplwai_aero_s2(klon), sollwai_aero_s2(klon))
     975      ALLOCATE(toplwad0_aero_s2(klon), sollwad0_aero_s2(klon))
     976
    921977! FH Ajout de celles necessaires au phys_output_write_mod
    922978
    … …  
    10381094      ALLOCATE(wfevap(klon, nbsrf))
    10391095      ALLOCATE(evap_pot(klon, nbsrf))
    1040 ! FC
     1096! FC 
    10411097      ALLOCATE(zxfluxq(klon,klev),zxfluxt(klon,klev))
    10421098!
    … …  
    10451101      ALLOCATE(pmflxr(klon, klev+1), pmflxs(klon, klev+1))
    10461102      ALLOCATE(wdtrainA(klon,klev),wdtrainS(klon,klev),wdtrainM(klon,klev))
    1047       ALLOCATE(dnwd(klon, klev), upwd(klon, klev) )
     1103      ALLOCATE(dnwd(klon, klev), upwd(klon, klev))
    10481104      ALLOCATE(ep(klon,klev))                          ! epmax_cape
    1049       ALLOCATE(da(klon,klev), mp(klon,klev) )
    1050       ALLOCATE(phi(klon,klev,klev) )
    1051       ALLOCATE(wght_cvfd(klon,klev) )
    1052       ALLOCATE(phi2(klon,klev,klev) )
     1105      ALLOCATE(da(klon,klev), mp(klon,klev))
     1106      ALLOCATE(phi(klon,klev,klev))
     1107      ALLOCATE(wght_cvfd(klon,klev))
     1108      ALLOCATE(phi2(klon,klev,klev))
    10531109      ALLOCATE(d1a(klon,klev), dam(klon,klev))
    1054       ALLOCATE(ev(klon,klev) )
    1055       ALLOCATE(elij(klon,klev,klev) )
    1056       ALLOCATE(qtaa(klon,klev) )
    1057       ALLOCATE(clw(klon,klev) )
    1058       ALLOCATE(epmlmMm(klon,klev,klev), eplaMm(klon,klev) )
    1059       ALLOCATE(sij(klon,klev,klev) )
     1110      ALLOCATE(ev(klon,klev))
     1111      ALLOCATE(elij(klon,klev,klev))
     1112      ALLOCATE(qtaa(klon,klev))
     1113      ALLOCATE(clw(klon,klev))
     1114      ALLOCATE(epmlmMm(klon,klev,klev), eplaMm(klon,klev))
     1115      ALLOCATE(sij(klon,klev,klev))
    10601116#ifdef ISO
    10611117      ALLOCATE(xtwdtrainA(ntraciso,klon,klev))
    … …  
    11281184
    11291185!--POPRECIP variables
    1130       ALLOCATE(dqreva(klon,klev),dqssub(klon,klev))
    1131       ALLOCATE(qraindiag(klon,klev), qsnowdiag(klon,klev))
     1186      ALLOCATE(qraindiag(klon,klev), qsnowdiag(klon,klev))
     1187      ALLOCATE(dqreva(klon,klev), dqssub(klon,klev))
    11321188      ALLOCATE(dqrauto(klon,klev), dqrcol(klon,klev), dqrmelt(klon,klev), dqrfreez(klon,klev))
    11331189      ALLOCATE(dqsauto(klon,klev), dqsagg(klon,klev), dqsrim(klon,klev), dqsmelt(klon,klev), dqsfreez(klon,klev))
    11341190
    11351191#ifdef CPP_StratAer
     1192      ALLOCATE (d_q_emiss(klon,klev))
    11361193      ALLOCATE (R2SO4(klon,klev))
    11371194      ALLOCATE (DENSO4(klon,klev))
    … …  
    11471204      ALLOCATE (OCS_lifetime(klon,klev))
    11481205      ALLOCATE (SO2_lifetime(klon,klev))
     1206      ALLOCATE (H2SO4_lifetime(klon,klev))
     1207      ALLOCATE (O3_clim(klon,klev))
    11491208      ALLOCATE (alpha_bin(nbands_sw_rrtm+nbands_lw_rrtm+nwave,nbtr))
    11501209      ALLOCATE (piz_bin(nbands_sw_rrtm+nbands_lw_rrtm+nwave,nbtr))
    … …  
    11801239USE indice_sol_mod
    11811240IMPLICIT NONE
    1182       DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri, qbs_seri)
     1241      DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri, qbs_seri,qx_seri)
    11831242      DEALLOCATE(u_seri,v_seri)
    11841243      DEALLOCATE(l_mixmin,l_mix,wprime)
    11851244      DEALLOCATE(pbl_eps)
    1186 
     1245      DEALLOCATE(rhcl)
    11871246      DEALLOCATE(tr_seri)
    11881247      DEALLOCATE(d_t_dyn,d_q_dyn)
    … …  
    12061265      DEALLOCATE(d_u_ajs,d_v_ajs)
    12071266      DEALLOCATE(d_t_eva,d_q_eva)
     1267      DEALLOCATE(d_qx_eva)
    12081268      DEALLOCATE(d_ql_eva,d_qi_eva)
    12091269      DEALLOCATE(d_t_lscst,d_q_lscst)
    … …  
    12141274      DEALLOCATE(d_t_bs,d_q_bs,d_qbs_bs)
    12151275#ifdef ISO
    1216       deallocate(xt_seri,xtl_seri,xts_seri)
     1276      deallocate(xt_seri,xtl_seri,xts_seri,xtbs_seri)
    12171277      DEALLOCATE(d_xtl_eva,d_xti_eva)
    1218       deallocate(d_xt_dyn,d_xtl_dyn,d_xts_dyn)
     1278      deallocate(d_xt_dyn,d_xtl_dyn,d_xts_dyn,d_xtbs_dyn)
    12191279      deallocate(d_xt_con)
    12201280      deallocate(d_xt_wake)
    … …  
    13061366      DEALLOCATE(solsw_aerop, solsw0_aerop)
    13071367      DEALLOCATE(topswcf_aerop, solswcf_aerop)
    1308 
     1368!AI Aerosols
     1369      DEALLOCATE(m_allaer)
    13091370!CK LW diagnostics
    13101371      DEALLOCATE(toplwad_aerop, sollwad_aerop)
    13111372      DEALLOCATE(toplwai_aerop, sollwai_aerop)
    13121373      DEALLOCATE(toplwad0_aerop, sollwad0_aerop)
     1374
     1375!AI Ajout pour Ecrad (3Deffect)
     1376      DEALLOCATE(topswad_aero_s2, solswad_aero_s2)
     1377      DEALLOCATE(topswai_aero_s2, solswai_aero_s2)
     1378      DEALLOCATE(topswad0_aero_s2, solswad0_aero_s2)
     1379      DEALLOCATE(topsw_aero_s2, topsw0_aero_s2)
     1380      DEALLOCATE(solsw_aero_s2, solsw0_aero_s2)
     1381      DEALLOCATE(topswcf_aero_s2, solswcf_aero_s2)
     1382!CK LW diagnostics
     1383      DEALLOCATE(toplwad_aero_s2, sollwad_aero_s2)
     1384      DEALLOCATE(toplwai_aero_s2, sollwai_aero_s2)
     1385      DEALLOCATE(toplwad0_aero_s2, sollwad0_aero_s2)     
    13131386
    13141387! FH Ajout de celles necessaires au phys_output_write_mod
    … …  
    13521425      DEALLOCATE(zxfqcalving, zxfluxlat)
    13531426      DEALLOCATE(zxrunofflic)
    1354       DEALLOCATE(zxustartlic, zxrhoslic)
     1427      DEALLOCATE(zxustartlic, zxrhoslic, zxqsaltlic)
    13551428      DEALLOCATE(zxtsol, snow_lsc, zxfqfonte, zxqsurf)
    13561429      DEALLOCATE(rain_lsc)
    … …  
    14201493      DEALLOCATE(upwd, dnwd)
    14211494      DEALLOCATE(ep)
    1422       DEALLOCATE(da, mp )
    1423       DEALLOCATE(phi )
    1424       DEALLOCATE(wght_cvfd )
    1425       DEALLOCATE(phi2 )
     1495      DEALLOCATE(da, mp)
     1496      DEALLOCATE(phi)
     1497      DEALLOCATE(wght_cvfd)
     1498      DEALLOCATE(phi2)
    14261499      DEALLOCATE(d1a, dam)
    1427       DEALLOCATE(ev )
    1428       DEALLOCATE(elij )
    1429       DEALLOCATE(qtaa )
    1430       DEALLOCATE(clw )
    1431       DEALLOCATE(epmlmMm, eplaMm )
    1432       DEALLOCATE(sij )
     1500      DEALLOCATE(ev)
     1501      DEALLOCATE(elij)
     1502      DEALLOCATE(qtaa)
     1503      DEALLOCATE(clw)
     1504      DEALLOCATE(epmlmMm, eplaMm)
     1505      DEALLOCATE(sij)
    14331506#ifdef ISO
    14341507      DEALLOCATE(xtwdtrainA)
    … …  
    14701543      DEALLOCATE(rneb)
    14711544      DEALLOCATE(pfraclr,pfracld)
     1545      DEALLOCATE (zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic)
    14721546      DEALLOCATE(distcltop)
    14731547      DEALLOCATE(temp_cltop)
    1474       DEALLOCATE (zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic)
    14751548#ifdef ISO
    14761549      DEALLOCATE (zxxtsnow,xtVprecip,xtVprecipi,pxtrfl,pxtsfl)
    … …  
    14931566
    14941567!--POPRECIP variables
    1495       DEALLOCATE(qraindiag, qsnowdiag) 
    1496       DEALLOCATE(dqreva,dqssub)
    1497       DEALLOCATE(dqrauto,dqrcol,dqrmelt,dqrfreez)
    1498       DEALLOCATE(dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez)
     1568      DEALLOCATE(qraindiag, qsnowdiag)
     1569      DEALLOCATE(dqreva, dqssub)
     1570      DEALLOCATE(dqrauto, dqrcol, dqrmelt, dqrfreez)
     1571      DEALLOCATE(dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
    14991572
    15001573#ifdef CPP_StratAer
    15011574! variables for strat. aerosol CK
     1575      DEALLOCATE (d_q_emiss)
    15021576      DEALLOCATE (R2SO4)
    15031577      DEALLOCATE (DENSO4)
    … …  
    15071581      DEALLOCATE (SO2_lifetime)
    15081582      DEALLOCATE (OCS_lifetime)
     1583      DEALLOCATE (H2SO4_lifetime)
     1584      DEALLOCATE (O3_clim)
    15091585      DEALLOCATE (alpha_bin)
    15101586      DEALLOCATE (piz_bin)

  • LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90

    r4976 r4982  
    3939    USE ioipsl_getin_p_mod, ONLY : getin_p
    4040    USE indice_sol_mod
    41     USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac
     41    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac,ivap,iliq,isol
    4242    USE readTracFiles_mod, ONLY: addPhase
    4343    USE strings_mod,  ONLY: strIdx
    … …  
    9494    USE phys_output_var_mod, ONLY : cloud_cover_sw, cloud_cover_sw_s2
    9595
     96
    9697    !USE cmp_seri_mod
    9798!    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
    … …  
    118119    USE chem_rep, ONLY: Init_chem_rep_xjour, d_q_rep, d_ql_rep, d_qi_rep, &
    119120                        ptrop, ttrop, ztrop, gravit, itroprep, Z1, Z2, fac, B
     121    USE strataer_local_var_mod
     122    USE strataer_emiss_mod, ONLY: strataer_emiss_init
    120123#endif
    121124#if defined INCA || defined REPROBUS
    … …  
    132135
    133136#ifdef CPP_StratAer
     137    USE phys_local_var_mod, ONLY: d_q_emiss
    134138    USE strataer_local_var_mod
    135139    USE strataer_nuc_mod, ONLY: strataer_nuc_init
    136140    USE strataer_emiss_mod, ONLY: strataer_emiss_init
    137141#endif
    138 
    139142
    140143    USE lmdz_xios, ONLY: xios_update_calendar, xios_context_finalize
    … …  
    154157        & modif_ratqs,essai_convergence,iso_init,ridicule_rain,tnat, &
    155158        & ridicule,ridicule_snow
    156     USE isotopes_routines_mod, ONLY: iso_tritium
     159    USE isotopes_routines_mod, ONLY: iso_tritium,dispatch,together
    157160#ifdef ISOVERIF
    158161    USE isotopes_verif_mod, ONLY: errmax,errmaxrel, &
    … …  
    189192!!!!!!!!!!!!!!!!!!  END "USE" for CPP keys !!!!!!!!!!!!!!!!!!!!!!
    190193
     194USE physiqex_mod, ONLY : physiqex
    191195USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
    192196       ! [Variables internes non sauvegardees de la physique]
    193197       ! Variables locales pour effectuer les appels en serie
    194198       t_seri,q_seri,ql_seri,qs_seri,qbs_seri,u_seri,v_seri,tr_seri,rneb_seri, &
     199       qx_seri, & ! CR
     200       rhcl, &       
    195201       ! Dynamic tendencies (diagnostics)
    196202       d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_qbs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn,d_rneb_dyn, &
    … …  
    207213       !
    208214       d_t_eva,d_q_eva,d_ql_eva,d_qi_eva, &
     215       d_qx_eva, &
    209216       d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc, &
    210217       d_t_lscst,d_q_lscst, &
    … …  
    245252       toplwai_aero,sollwai_aero,   &
    246253       toplwad0_aero,sollwad0_aero, &
     254       !pour Ecrad
     255       topswad_aero_s2, solswad_aero_s2,   &
     256       topswai_aero_s2, solswai_aero_s2,   &
     257       topswad0_aero_s2, solswad0_aero_s2, &
     258       topsw_aero_s2, topsw0_aero_s2,      &
     259       solsw_aero_s2, solsw0_aero_s2,      &
     260       topswcf_aero_s2, solswcf_aero_s2,   &
     261       !LW diagnostics
     262       toplwad_aero_s2, sollwad_aero_s2,   &
     263       toplwai_aero_s2, sollwai_aero_s2,   &
     264       toplwad0_aero_s2, sollwad0_aero_s2, &
    247265       !
    248266       topsw_aero,solsw_aero,       &
    … …  
    263281       toplwai_aerop, sollwai_aerop,   &
    264282       toplwad0_aerop, sollwad0_aerop, &
     283       !pour Ecrad
     284       topswad_aero_s2, solswad_aero_s2,   &
     285       topswai_aero_s2, solswai_aero_s2,   &
     286       topswad0_aero_s2, solswad0_aero_s2, &
     287       topsw_aero_s2, topsw0_aero_s2,      &
     288       solsw_aero_s2, solsw0_aero_s2,      &
     289       topswcf_aero_s2, solswcf_aero_s2,   &
     290       !LW diagnostics
     291       toplwad_aero_s2, sollwad_aero_s2,   &
     292       toplwai_aero_s2, sollwai_aero_s2,   &
     293       toplwad0_aero_s2, sollwad0_aero_s2, &
    265294       !
    266295       ptstar, pt0, slp, &
    … …  
    271300       JrNt,                             &
    272301       dthmin, evap, snowerosion,fder, plcl, plfc,   &
    273        prw, prlw, prsw, prbsw,                  &
     302       prw, prlw, prsw, prbsw, water_budget,         &
    274303       s_lcl, s_pblh, s_pblt, s_therm,   &
    275304       cdragm, cdragh,                   &
    … …  
    359388       t2m, fluxlat,  &
    360389       fsollw, evap_pot,  &
    361        fsolsw, wfbils, wfevap,  &
     390       fsolsw, wfbils, wfevap, & 
    362391       prfl, psfl,bsfl, fraca, Vprecip,  &
    363392       zw2,  &
    … …  
    373402       rneb,  &
    374403       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic, &
    375        zxfluxt,zxfluxq
    376 
    377 
    378 #ifdef ISO
    379        USE phys_local_var_mod, ONLY: xt_seri,xtl_seri,xts_seri, &
     404       zxfluxt,zxfluxq 
     405
     406
     407#ifdef ISO
     408       USE phys_local_var_mod, ONLY: xt_seri,xtl_seri,xts_seri,xtbs_seri, &
    380409       d_xt_eva,d_xtl_eva,d_xti_eva,           &
    381        d_xt_dyn,d_xtl_dyn,d_xts_dyn,          &
     410       d_xt_dyn,d_xtl_dyn,d_xts_dyn,d_xtbs_dyn, &
    382411       d_xt_con, d_xt_wake,                    &
    383412       d_xt_ajsb, d_xt_ajs,                    &
    … …  
    405434       USE phys_output_var_mod, ONLY: scdnc, cldncl, reffclwtop, lcc, reffclws, &
    406435       reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra
     436       USE output_physiqex_mod, ONLY: output_physiqex
    407437
    408438
    … …  
    549579    !
    550580    ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional), blowing snow (optional)
    551     INTEGER,SAVE :: ivap, iliq, isol, irneb, ibs
    552 !$OMP THREADPRIVATE(ivap, iliq, isol, irneb, ibs)
    553     !
     581!    INTEGER,SAVE :: ivap, iliq, isol, irneb, ibs
     582!!$OMP THREADPRIVATE(ivap, iliq, isol, irneb, ibs)
     583! Camille Risi 25 juillet 2023: ivap,iliq,isol deja definis dans infotrac_phy.
     584! Et ils sont utiles ailleurs que dans physiq_mod (ex:
     585! reevap -> je commente les 2 lignes au dessus et je laisse la definition
     586! plutot dans infotrac_phy
     587    INTEGER,SAVE :: irneb, ibs
     588!$OMP THREADPRIVATE(irneb, ibs)
     589!
    554590    !
    555591    ! Variables argument:
    … …  
    805841    real therm_tke_max0(klon)   ! TKE dans les thermiques au LCL
    806842    real env_tke_max0(klon)     ! TKE dans l'environnement au LCL
     843    INTEGER, SAVE :: iflag_thermcell_tke ! transtport TKE by thermals
     844    !$OMP THREADPRIVATE(iflag_thermcell_tke)
    807845
    808846!JLD    !---D\'eclenchement stochastique
    … …  
    893931    EXTERNAL ajsec     ! ajustement sec
    894932    EXTERNAL conlmd    ! convection (schema LMD)
    895     !KE43
    896933    EXTERNAL conema3  ! convect4.3
    897     !AA
    898     ! JBM (3/14) fisrtilp_tr not loaded
    899     ! EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
    900     !                          ! stockage des coefficients necessaires au
    901     !                          ! lessivage OFF-LINE et ON-LINE
    902934    EXTERNAL hgardfou  ! verifier les temperatures
    903935    EXTERNAL nuage     ! calculer les proprietes radiatives
    … …  
    921953    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    922954    !
    923     REAL rhcl(klon,klev)    ! humiditi relative ciel clair
     955!    REAL rhcl(klon,klev)    ! humiditi relative ciel clair
    924956    REAL dialiq(klon,klev)  ! eau liquide nuageuse
    925957    REAL diafra(klon,klev)  ! fraction nuageuse
    … …  
    953985    REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
    954986    !
    955 #ifdef INCA
    956     REAL zxsnow_dummy(klon)
    957 #endif
    958987    REAL zsav_tsol(klon)
    959988    !
    … …  
    10611090    !$OMP THREADPRIVATE(ok_bug_split_th)
    10621091
     1092    ! Logical switch to a bug : modifying directly wake_deltat  by adding
     1093    ! the (w) dry adjustment tendency to wake_deltat
     1094    LOGICAL, SAVE :: ok_bug_ajs_cv = .TRUE.
     1095    !$OMP THREADPRIVATE(ok_bug_ajs_cv)
    10631096    !
    10641097    !********************************************************
    … …  
    12651298    ! Declarations pour Simulateur COSP
    12661299    !============================================================
     1300    ! AI 10-22
     1301#ifdef CPP_COSP
     1302    include "ini_COSP.h"
     1303#endif
     1304#ifdef CPP_COSPV2
     1305    include "ini_COSP.h"
     1306#endif
    12671307    real :: mr_ozone(klon,klev), phicosp(klon,klev)
    12681308
    … …  
    13401380
    13411381    REAL, dimension(klon,klev) :: t_env,q_env
     1382#ifdef ISO
     1383    real, dimension(ntraciso,klon,klev) ::  xt_env
     1384#endif
    13421385
    13431386    REAL, dimension(klon) :: pr_et
    … …  
    13501393    !AI namelist pour gerer le double appel de Ecrad
    13511394    CHARACTER(len=512) :: namelist_ecrad_file
     1395
     1396    !======================================================================!
     1397    ! Bifurcation vers un nouveau moniteur physique pour experimenter      !
     1398    ! des solutions et préparer le couplage avec la physique de MesoNH     !
     1399    ! 14 mai 2023                                                          !
     1400    !======================================================================!
     1401    if (debut) then                                                        !
     1402       iflag_physiq=0
     1403       call getin_p('iflag_physiq', iflag_physiq)                          !
     1404    endif                                                                  !
     1405    if ( iflag_physiq == 2 ) then
     1406#ifdef ISO
     1407       abort_message='isotopes pas encore dans physiqex'
     1408       CALL abort_physic(modname,abort_message,1)
     1409#endif
     1410       call physiqex (nlon,nlev, &                                         !
     1411       debut,lafin,pdtphys_, &                                             !
     1412       paprs,pplay,pphi,pphis,presnivs, &                                  !
     1413       u,v,rot,t,qx, &                                                     !
     1414       flxmass_w, &                                                        !
     1415       d_u, d_v, d_t, d_qx, d_ps)                                          !
     1416       return                                                              !
     1417    endif                                                                  !
     1418    !======================================================================!
     1419
    13521420
    13531421    pi = 4. * ATAN(1.)
    … …  
    13661434    phys_tstep=NINT(pdtphys)
    13671435    IF (.NOT. using_xios) missing_val=nf90_fill_real
    1368 #ifdef CPP_XIOS
    1369 ! switch to XIOS LMDZ physics context
    1370     IF (.NOT. debut .AND. is_omp_master) THEN
    1371        CALL wxios_set_context()
    1372        CALL xios_update_calendar(itap+1)
     1436
     1437    IF (using_xios) THEN
     1438      ! switch to XIOS LMDZ physics context
     1439      IF (.NOT. debut .AND. is_omp_master) THEN
     1440        CALL wxios_set_context()
     1441        CALL xios_update_calendar(itap+1)
     1442      ENDIF
    13731443    ENDIF
    1374 #endif
    13751444
    13761445    !======================================================================
    … …  
    13781447    ! Utilise notamment en 1D mais peut etre active egalement en 3D
    13791448    ! en imposant la valeur de igout.
    1380     !======================================================================d
     1449    !======================================================================
    13811450    IF (prt_level.ge.1) THEN
    13821451       igout=klon/2+1/klon
    … …  
    14101479       irneb= strIdx(tracers(:)%name, addPhase('H2O', 'r'))
    14111480       ibs  = strIdx(tracers(:)%name, addPhase('H2O', 'b'))
    1412        CALL init_etat0_limit_unstruct
    1413        IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
     1481!       CALL init_etat0_limit_unstruct
     1482!       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
    14141483       !CR:nvelles variables convection/poches froides
    14151484
    … …  
    14341503            read_climoz, &
    14351504            alp_offset)
     1505       CALL init_etat0_limit_unstruct
     1506       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
    14361507       CALL phys_state_var_init(read_climoz)
    14371508       CALL phys_output_var_init
    14381509       IF (read_climoz>=1 .AND. create_etat0_limit .AND. grid_type==unstructured) &
    14391510          CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
     1511
     1512#ifdef REPROBUS
     1513       CALL strataer_init
     1514       CALL strataer_emiss_init
     1515#endif
    14401516
    14411517#ifdef CPP_StratAer
    … …  
    14811557
    14821558        IF (ok_bs) THEN
     1559#ifdef ISO
    14831560          abort_message='blowing snow cannot be activated with water isotopes yet'
    14841561          CALL abort_physic(modname,abort_message, 1)
     1562#endif
    14851563         IF ((ok_ice_sursat.AND.nqo .LT.5).OR.(.NOT.ok_ice_sursat.AND.nqo.LT.4)) THEN
    14861564             WRITE (lunout, *) 'activation of blowing snow needs a specific H2O tracer', &
    … …  
    14901568         ENDIF
    14911569        ENDIF
     1570
    14921571       Ncvpaseq1 = 0
    14931572       dnwd0=0.0
    … …  
    15311610       tau_gl=86400.*tau_gl
    15321611       WRITE(lunout,*) 'debut physiq_mod tau_gl=',tau_gl
     1612       iflag_thermcell_tke=0
     1613       call getin_p('iflag_thermcell_tke', iflag_thermcell_tke)                          !
    15331614
    15341615       CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond)
    … …  
    15531634       CALL getin_p('ok_bug_cv_trac',ok_bug_cv_trac)
    15541635       CALL getin_p('ok_bug_split_th',ok_bug_split_th)
     1636       CALL getin_p('ok_bug_ajs_cv',ok_bug_ajs_cv)
    15551637       fl_ebil = 0 ! by default, conservation diagnostics are desactivated
    15561638       CALL getin_p('fl_ebil',fl_ebil)
    … …  
    15891671       CALL infocfields_init
    15901672
     1673       !AI 08 2023
    15911674#ifdef CPP_ECRAD
    15921675       ok_3Deffect=.false.
    … …  
    18681951      IF (.NOT. create_etat0_limit) CALL init_readaerosolstrato(flag_aerosol_strat)  !! initialise aero strato from file for XIOS interpolation (unstructured_grid)
    18691952
     1953      if (ok_cosp) then
    18701954#ifdef CPP_COSP
    1871       IF (ok_cosp) THEN
    1872 !           DO k = 1, klev
    1873 !             DO i = 1, klon
    1874 !               phicosp(i,k) = pphi(i,k) + pphis(i)
    1875 !             ENDDO
    1876 !           ENDDO
     1955        ! A.I : Initialisations pour le 1er passage a Cosp
     1956        CALL ini_COSP(ref_liq_cosp0,ref_ice_cosp0,pctsrf_cosp0,zu10m_cosp0,zv10m_cosp0, &
     1957               zxtsol_cosp0,zx_rh_cosp0,cldfra_cosp0,rnebcon_cosp0,flwc_cosp0, &
     1958               fiwc_cosp0,prfl_cosp0,psfl_cosp0,pmflxr_cosp0,pmflxs_cosp0, &
     1959               mr_ozone_cosp0,cldtau_cosp0,cldemi_cosp0,JrNt_cosp0)
     1960
    18771961        CALL phys_cosp(itap,phys_tstep,freq_cosp, &
     1962               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
     1963               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
     1964               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
     1965               JrNt_cosp0,ref_liq_cosp0,ref_ice_cosp0, &
     1966               pctsrf_cosp0, &
     1967               zu10m_cosp0,zv10m_cosp0,pphis, &
     1968               pphi,paprs(:,1:klev),pplay,zxtsol_cosp0,t, &
     1969               qx(:,:,ivap),zx_rh_cosp0,cldfra_cosp0,rnebcon_cosp0,flwc_cosp0,fiwc_cosp0, &
     1970               prfl_cosp0(:,1:klev),psfl_cosp0(:,1:klev), &
     1971               pmflxr_cosp0(:,1:klev),pmflxs_cosp0(:,1:klev), &
     1972               mr_ozone_cosp0,cldtau_cosp0, cldemi_cosp0)
     1973#endif
     1974
     1975#ifdef CPP_COSP2
     1976        CALL ini_COSP(ref_liq_cosp0,ref_ice_cosp0,pctsrf_cosp0,zu10m_cosp0,zv10m_cosp0, &
     1977               zxtsol_cosp0,zx_rh_cosp0,cldfra_cosp0,rnebcon_cosp0,flwc_cosp0, &
     1978               fiwc_cosp0,prfl_cosp0,psfl_cosp0,pmflxr_cosp0,pmflxs_cosp0, &
     1979               mr_ozone_cosp0,cldtau_cosp0,cldemi_cosp0,JrNt_cosp0)
     1980     
     1981        CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
    18781982               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
    18791983               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
    … …  
    18871991               pmflxr(:,1:klev),pmflxs(:,1:klev), &
    18881992               mr_ozone,cldtau, cldemi)
    1889       ENDIF
    1890 #endif
    1891 
    1892 #ifdef CPP_COSP2
    1893         IF (ok_cosp) THEN
    1894 !           DO k = 1, klev
    1895 !             DO i = 1, klon
    1896 !               phicosp(i,k) = pphi(i,k) + pphis(i)
    1897 !             ENDDO
    1898 !           ENDDO
    1899           CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
    1900                ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
    1901                ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
    1902                klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
    1903                JrNt,ref_liq,ref_ice, &
    1904                pctsrf(:,is_ter)+pctsrf(:,is_lic), &
    1905                zu10m,zv10m,pphis, &
    1906                zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
    1907                qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
    1908                prfl(:,1:klev),psfl(:,1:klev), &
    1909                pmflxr(:,1:klev),pmflxs(:,1:klev), &
    1910                mr_ozone,cldtau, cldemi)
    1911        ENDIF
    19121993#endif
    19131994
    19141995#ifdef CPP_COSPV2
    1915         IF (ok_cosp) THEN
    1916            DO k = 1, klev
    1917              DO i = 1, klon
    1918                phicosp(i,k) = pphi(i,k) + pphis(i)
    1919              ENDDO
    1920            ENDDO
    19211996          CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
    19221997               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
    … …  
    19312006               pmflxr(:,1:klev),pmflxs(:,1:klev), &
    19322007               mr_ozone,cldtau, cldemi)
    1933        ENDIF
    1934 #endif
     2008#endif
     2009      ENDIF
    19352010
    19362011       !
    … …  
    20062081
    20072082
    2008 #ifdef CPP_XIOS
    2009        IF (is_omp_master) CALL xios_update_calendar(1)
    2010 #endif
     2083       IF (using_xios) THEN
     2084         IF (is_omp_master) CALL xios_update_calendar(1)
     2085       ENDIF
     2086       
    20112087       IF(read_climoz>=1 .AND. create_etat0_limit) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
    20122088       CALL create_etat0_limit_unstruct
    … …  
    22112287        IF (is_omp_master) CALL xios_get_field_attr("temp",default_value=missing_val)
    22122288        CALL bcast_omp(missing_val)
    2213        
     2289
    22142290       !
    22152291       ! Now we activate some double radiation call flags only if some
    … …  
    25392615          u_seri(i,k)  = u(i,k)
    25402616          v_seri(i,k)  = v(i,k)
     2617          qx_seri(i,k,:)  = qx(i,k,:)
    25412618          q_seri(i,k)  = qx(i,k,ivap)
    25422619          ql_seri(i,k) = qx(i,k,iliq)
    … …  
    25522629             qs_seri(i,k) = qx(i,k,isol)
    25532630             IF (ok_ice_sursat) THEN
    2554              rneb_seri(i,k) = qx(i,k,irneb)
     2631               rneb_seri(i,k) = qx(i,k,irneb)
    25552632             ENDIF
    25562633             IF (ok_bs) THEN
    2557              qbs_seri(i,k)= qx(i,k,ibs)
     2634               qbs_seri(i,k)= qx(i,k,ibs)
    25582635             ENDIF
    2559 
    25602636          ENDIF
    2561 
    2562 
    25632637       ENDDO
    25642638    ENDDO
    … …  
    25822656      DO k = 1, klev
    25832657       DO i = 1, klon
     2658          xtbs_seri(ixt,i,k) = 0.
    25842659          xt_seri(ixt,i,k)  = qx(i,k,iqIsoPha(ixt,ivap))
    25852660          xtl_seri(ixt,i,k) = qx(i,k,iqIsoPha(ixt,iliq))
    … …  
    26022677    qql1(:)=0.0
    26032678    DO k = 1, klev
    2604       qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k)+qbs_seri(:,k))*zmasse(:,k)
     2679      qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k))*zmasse(:,k)
     2680      IF (nqo >= 3) THEN
     2681        qql1(:)=qql1(:)+qs_seri(:,k)*zmasse(:,k)
     2682      ENDIF
     2683      IF (ok_bs) THEN
     2684        qql1(:)=qql1(:)+qbs_seri(:,k)*zmasse(:,k)
     2685      ENDIF
    26052686    ENDDO
    26062687#ifdef ISO
    2607 #ifdef ISOVERIF
    2608         write(*,*) 'physiq tmp 1913'
    2609 #endif
    26102688    do ixt=1,ntraciso
    26112689    xtql1(ixt,:)=0.0
    26122690    DO k = 1, klev
    2613       xtql1(ixt,:)=xtql1(ixt,:)+(xt_seri(ixt,:,k)+xtl_seri(ixt,:,k)+xts_seri(ixt,:,k))*zmasse(:,k)
    2614     ENDDO
     2691      xtql1(ixt,:)=xtql1(ixt,:)+(xt_seri(ixt,:,k)+xtl_seri(ixt,:,k))*zmasse(:,k)
     2692      IF (nqo >= 3) THEN
     2693        xtql1(ixt,:)=xtql1(ixt,:)+xts_seri(ixt,:,k)*zmasse(:,k)
     2694      ENDIF
     2695      IF (ok_bs) THEN
     2696        xtql1(ixt,:)=xtql1(ixt,:)+xtbs_seri(ixt,:,k)*zmasse(:,k)
     2697      ENDIF
     2698    ENDDO !DO k = 1, klev
    26152699    enddo !do ixt=1,ntraciso
    26162700#endif
    … …  
    26252709         IF(.NOT.tracers(iq)%isInPhysics) CYCLE
    26262710         itr = itr+1
    2627 !#ifdef ISOVERIF
    2628 !         write(*,*) 'physiq 1973: itr,iq=',itr,iq
    2629 !         write(*,*) 'qx(1,1,iq)=',qx(1,1,iq)
    2630 !#endif
    2631          DO  k = 1, klev
     2711          DO  k = 1, klev
    26322712             DO  i = 1, klon
    26332713                tr_seri(i,k,itr) = qx(i,k,iq)
    … …  
    27442824              d_xts_dyn(ixt,i,k) =  &
    27452825     &           (xts_seri(ixt,i,k)-xts_ancien(ixt,i,k))/phys_tstep
     2826              d_xtbs_dyn(ixt,i,k) =  &
     2827     &           (xtbs_seri(ixt,i,k)-xtbs_ancien(ixt,i,k))/phys_tstep
    27462828            enddo ! do ixt=1,ntraciso
    27472829         ENDDO
    … …  
    27572839           call iso_verif_noNaN(d_xtl_dyn(ixt,i,k),'physiq 2220d')
    27582840           call iso_verif_noNaN(d_xts_dyn(ixt,i,k),'physiq 2220e')
     2841           call iso_verif_noNaN(d_xtbs_dyn(ixt,i,k),'physiq 2220f')
    27592842           enddo ! do ixt=1,ntraciso
    27602843         enddo
    … …  
    28392922       d_qbs_dyn(:,:)=0.0
    28402923       ancien_ok = .TRUE.
     2924#ifdef ISO
     2925       d_xtbs_dyn(:,:,:)=0.0
     2926#endif
    28412927    ENDIF
    28422928    !
    … …  
    29623048    ! Re-evaporer l'eau liquide nuageuse
    29633049    !
    2964      CALL reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
    2965    &         d_t_eva,d_q_eva,d_ql_eva,d_qi_eva &
    2966 #ifdef ISO
    2967              ,xt_seri,xtl_seri,xts_seri,d_xt_eva,d_xtl_eva,d_xti_eva &
    2968 #endif
    2969    &     )
     3050     CALL reevap (klon,klev,iflag_ice_thermo,t_seri,qx_seri, &
     3051   &         d_t_eva,d_qx_eva)
     3052
     3053     call dispatch(klon,klev,qx_seri,q_seri,xt_seri,ql_seri,xtl_seri,qs_seri,xts_seri)
     3054     call dispatch(klon,klev,d_qx_eva,d_q_eva,d_xt_eva,d_ql_eva,d_xtl_eva,d_qi_eva,d_xti_eva)
     3055
     3056
     3057#ifdef ISO
     3058#ifdef ISOVERIF
     3059 DO k = 1, klev
     3060     DO i = 1, klon
     3061      do ixt=1,ntraciso
     3062      call iso_verif_noNaN(xt_seri(ixt,i,k), &
     3063     &     'reevap 2417: apres evap tot')
     3064      enddo
     3065      if (iso_eau.gt.0) then
     3066              call iso_verif_egalite_choix( &
     3067     &           xt_seri(iso_eau,i,k),q_seri(i,k), &
     3068     &          'reevap 1891, après réévap totale',errmax,errmaxrel)
     3069              call iso_verif_egalite_choix( &
     3070     &           xtl_seri(iso_eau,i,k),ql_seri(i,k), &
     3071     &          'reevap 2209, après réévap totale',errmax,errmaxrel)
     3072              call iso_verif_egalite_choix( &
     3073     &           xts_seri(iso_eau,i,k),qs_seri(i,k), &
     3074     &          'reevap 2209b, après réévap totale',errmax,errmaxrel)
     3075       endif !if (iso_eau.gt.0) then
     3076     
     3077      if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then 
     3078            if (q_seri(i,k).gt.ridicule) then 
     3079               if (iso_verif_o18_aberrant_nostop( &
     3080     &              xt_seri(iso_HDO,i,k)/q_seri(i,k), &
     3081     &              xt_seri(iso_O18,i,k)/q_seri(i,k), &
     3082     &              'reevap 2408: apres reevap totale').eq.1) then
     3083                  write(*,*) 'i,k,q_seri(i,k)=',i,k,q_seri(i,k)                       
     3084                  stop
     3085              endif !  if (iso_verif_o18_aberrant_nostop
     3086            endif !if (q_seri(i,k).gt.errmax) then   
     3087        endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then       
     3088#ifdef ISOTRAC     
     3089             call iso_verif_traceur(xt_seri(1,i,k), &
     3090     &           'reevap 2165c')
     3091             call iso_verif_traceur_pbidouille(xt_seri(1,i,k), &
     3092     &           'reevap 2165d')
     3093#endif
     3094       ENDDO
     3095    ENDDO               
     3096#endif                 
     3097#endif
     3098
    29703099
    29713100     CALL add_phys_tend &
    … …  
    30993228    ! Calcul de l'humidite de saturation au niveau du sol
    31003229
     3230! Tests Fredho, instensibilite au pas de temps -------------------------------
     3231! A detruire en 2024 une fois les tests documentes et les choix faits        !
     3232! Conservation des variables avant l'appel à l a diffusion pour les tehrmic  !
     3233    if (iflag_thermals_tenv / 10 == 1 ) then                                 !
     3234        do k=1,klev                                                          !
     3235           do i=1,klon                                                       !
     3236              t_env(i,k)=t_seri(i,k)                                         !
     3237              q_env(i,k)=q_seri(i,k)   
     3238#ifdef ISO
     3239              do ixt=1,ntraciso
     3240                xt_env(ixt,i,k)=xt_seri(ixt,i,k)
     3241              enddo
     3242#endif
     3243           enddo                                                             !
     3244        enddo                                                                !
     3245    else if (iflag_thermals_tenv / 10 == 2 ) then                            !
     3246        do k=1,klev                                                          !
     3247           do i=1,klon                                                       !
     3248              t_env(i,k)=t_seri(i,k)                                         !
     3249           enddo                                                             !
     3250        enddo                                                                !
     3251    endif                                                                    !
     3252! Tests Fredho, instensibilite au pas de temps -------------------------------
    31013253
    31023254
    … …  
    32873439          d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
    32883440          CALL add_wake_tend &
    3289              (d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy &
     3441             (d_deltat_vdf, d_deltaq_vdf, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy &
    32903442#ifdef ISO
    32913443               ,d_deltaxt_vdf &
    … …  
    33203472     &   )
    33213473       ENDIF
    3322 #ifdef ISOVERIF
    3323         write(*,*) 'physiq tmp 2736'
    3324 #endif
    3325 
    33263474       CALL prt_enerbil('vdf',itap)
     3475
    33273476       !--------------------------------------------------------------------
    33283477
    … …  
    36893838                ENDDO
    36903839             ENDDO
    3691              IF (iflag_adjwk == 2) THEN
     3840             IF (iflag_adjwk == 2 .AND. OK_bug_ajs_cv) THEN
    36923841               CALL add_wake_tend &
    3693                  (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy &
     3842                 (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy &
    36943843#ifdef ISO
    36953844                      ,d_deltaxt_ajs_cv  &
    36963845#endif
    36973846                 )
    3698              ENDIF  ! (iflag_adjwk == 2)
     3847             ENDIF  ! (iflag_adjwk == 2 .AND. OK_bug_ajs_cv)
    36993848          ENDIF  ! (iflag_adjwk >= 1)
    37003849       ENDIF ! (iflag_wake>=1)
    … …  
    44004549       !  ====
    44014550       IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
     4551       WRITE(lunout,*) 'WARNING : running without dry convection. Somme intermediate variables are not properly defined in physiq_mod.F90'
     4552       ! Reprendre proprement les initialisation ci dessouds si on veut vraiment utiliser l'option (FH)
     4553          fraca(:,:)=0.
     4554          fm_therm(:,:)=0.
     4555          ztv(:,:)=t_seri(:,:)
     4556          zqasc(:,:)=q_seri(:,:)
     4557          ztla(:,:)=0.
     4558          zthl(:,:)=0.
     4559          zpspsk(:,:)=(pplay(:,:)/100000.)**RKAPPA
    44024560
    44034561
    … …  
    44914649
    44924650       IF (iflag_thermals>=1) THEN
     4651
     4652! Tests Fredho, instensibilite au pas de temps -------------------------------
     4653! A detruire en 2024 une fois les tests documentes et les choix faits        !
     4654          if (iflag_thermals_tenv /10 == 0 ) then                            !
     4655            do k=1,klev                                                      !
     4656               do i=1,klon                                                   !
     4657                  t_env(i,k)=t_seri(i,k)                                     !
     4658                  q_env(i,k)=q_seri(i,k)                                     !
     4659#ifdef ISO
     4660                  do ixt=1,ntraciso
     4661                        xt_env(ixt,i,k)=xt_seri(ixt,i,k)
     4662                  enddo
     4663#endif
     4664               enddo                                                         !
     4665            enddo                                                            !
     4666          else if (iflag_thermals_tenv / 10 == 2 ) then                      !
     4667            do k=1,klev                                                      !
     4668               do i=1,klon                                                   !
     4669                  q_env(i,k)=q_seri(i,k)                                     !
     4670#ifdef ISO
     4671                  do ixt=1,ntraciso
     4672                        xt_env(ixt,i,k)=xt_seri(ixt,i,k)
     4673                  enddo
     4674#endif
     4675               enddo                                                         !
     4676            enddo                                                            !
     4677          else if (iflag_thermals_tenv / 10 == 3 ) then                      !
     4678            do k=1,klev                                                      !
     4679               do i=1,klon                                                   !
     4680                  t_env(i,k)=t(i,k)                                          !
     4681                  q_env(i,k)=qx(i,k,1)                                       !
     4682#ifdef ISO
     4683                  do ixt=1,ntraciso
     4684                        xt_env(ixt,i,k)=xt_seri(ixt,i,k)
     4685                  enddo
     4686#endif
     4687               enddo                                                         !
     4688            enddo                                                            !
     4689          endif                                                              !
     4690! Tests Fredho, instensibilite au pas de temps ------------------------------
     4691
    44934692          !jyg<
    44944693!!       IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
    … …  
    44994698                   t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
    45004699                   q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
     4700                   t_env(i,k)   = t_env(i,k) - wake_s(i)*wake_deltat(i,k)
     4701                   q_env(i,k)   = q_env(i,k) - wake_s(i)*wake_deltaq(i,k)
    45014702                   u_therm(i,k) = u_seri(i,k)
    45024703                   v_therm(i,k) = v_seri(i,k)
    … …  
    45044705                   do ixt=1,ntraciso
    45054706                     xt_therm(ixt,i,k) = xt_seri(ixt,i,k) - wake_s(i)*wake_deltaxt(ixt,i,k)
     4707                     xt_env(ixt,i,k) = xt_env(ixt,i,k) - wake_s(i)*wake_deltaxt(ixt,i,k)
    45064708                   enddo !do ixt=1,ntraciso
    45074709#endif
    … …  
    45464748               ,pplay,paprs,pphi,weak_inversion &
    45474749                        ! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
    4548                ,u_therm,v_therm,t_therm,q_therm,t_therm,q_therm,zqsat,debut &  !jyg
     4750               ,u_therm,v_therm,t_therm,q_therm,t_env,q_env,zqsat,debut &  !jyg
    45494751               ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
    45504752               ,fm_therm,entr_therm,detr_therm &
    … …  
    45654767               ,zqla,ztva &
    45664768#ifdef ISO         
    4567      &      ,xt_therm,d_xt_ajs &
     4769     &      ,xt_env,d_xt_ajs &
    45684770#ifdef DIAGISO         
    45694771     &      ,q_the,xt_the &
    … …  
    46004802             IF (ok_bug_split_th) THEN
    46014803               CALL add_wake_tend &
    4602                    (d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy &
     4804                   (d_deltat_the, d_deltaq_the, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy &
    46034805#ifdef ISO
    46044806                   ,d_deltaxt_the &
    … …  
    46074809             ELSE
    46084810               CALL add_wake_tend &
    4609                    (d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wake_k, 'the', abortphy &
     4811                   (d_deltat_the, d_deltaq_the, dsig0, dsig0, ddens0, ddens0, wake_k, 'the', abortphy &
    46104812#ifdef ISO
    46114813                   ,d_deltaxt_the &
    … …  
    46364838          ! Transport de la TKE par les panaches thermiques.
    46374839          ! FH : 2010/02/01
    4638           !     if (iflag_pbl.eq.10) then
    4639           !     call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
    4640           !    s           rg,paprs,pbl_tke)
    4641           !     endif
     4840               if (iflag_thermcell_tke==1) then
     4841               call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,rg,paprs,pbl_tke)
     4842               endif
    46424843          ! -------------------------------------------------------------------
    46434844
    … …  
    48955096    ELSE
    48965097
     5098    ! Camille Risi mai 2024: on ne met pas à jour ici pour ne pas s'mbêter à modifier fisrtilp
    48975099    CALL fisrtilp(phys_tstep,paprs,pplay, &
    48985100         t_seri, q_seri,ptconv,ratqs, &
    … …  
    55605762                !
    55615763             ENDIF
     5764<<<<<<< .mine
     5765          ELSE IF (iflag_rrtm .EQ.2) THEN    ! ecrad RADIATION
     5766#ifdef CPP_ECRAD
     5767             !--climatologies or INCA aerosols
     5768             CALL readaerosol_optic_ecrad( debut, aerosol_couple, ok_alw, ok_volcan, &
     5769                  flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
     5770                  pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
     5771                  tr_seri, mass_solu_aero, mass_solu_aero_pi) 
     5772#else
     5773                abort_message='You should compile with -rad ecrad if running with iflag_rrtm=2'
     5774                CALL abort_physic(modname,abort_message,1)
     5775#endif
     5776||||||| .r4942
     5777=======
    55625778          ELSE IF (iflag_rrtm .EQ.2) THEN    ! ecrad RADIATION
    55635779#ifdef CPP_ECRAD
    … …  
    55715787                CALL abort_physic(modname,abort_message,1)
    55725788#endif
     5789>>>>>>> .r4981
    55735790          ENDIF
     5791
    55745792       ELSE   !--flag_aerosol = 0
    55755793          tausum_aero(:,:,:) = 0.
    … …  
    58736091                ENDIF
    58746092                !
    5875                 ! AI namelist utilise pour l appel principal de radlwsw (ecrad)
    58766093                namelist_ecrad_file='namelist_ecrad'
    58776094                !
    … …  
    59176134                     cloud_cover_sw)
    59186135          ENDIF !ok_4xCO2atm
     6136
     6137! A.I aout 2023
     6138! Effet 3D des nuages Ecrad
     6139! a passer : nom du ficher namelist et cles ok_3Deffect
     6140! a declarer comme iflag_rrtm et a lire dans physiq.def
     6141#ifdef CPP_ECRAD
     6142          IF (ok_3Deffect) then
     6143!                print*,'ok_3Deffect = ',ok_3Deffect 
     6144                namelist_ecrad_file='namelist_ecrad_s2'
     6145                CALL radlwsw &
     6146                     (debut, dist, rmu0, fract,  &
     6147                     paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
     6148                     t_seri,q_seri,wo, &
     6149                     cldfrarad, cldemirad, cldtaurad, &
     6150                     ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie,  ok_volcan, flag_volc_surfstrat, &
     6151                     flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
     6152                     tau_aero, piz_aero, cg_aero, &
     6153                     tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
     6154                     tau_aero_lw_rrtm, &
     6155                     cldtaupi, &
     6156                     zqsat, flwc, fiwc, &
     6157                     ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
     6158                     namelist_ecrad_file, &
     6159! A modifier             
     6160                     heat_s2,heat0_s2,cool_s2,cool0_s2,albpla_s2, &
     6161                     heat_volc,cool_volc, &
     6162                     topsw_s2,toplw_s2,solsw_s2,solswfdiff_s2,sollw_s2, &
     6163                     sollwdown_s2, &
     6164                     topsw0_s2,toplw0_s2,solsw0_s2,sollw0_s2, &
     6165                     lwdnc0_s2, lwdn0_s2, lwdn_s2, lwupc0_s2, lwup0_s2, lwup_s2,  &
     6166                     swdnc0_s2, swdn0_s2, swdn_s2, swupc0_s2, swup0_s2, swup_s2, &
     6167                     topswad_aero_s2, solswad_aero_s2, &
     6168                     topswai_aero_s2, solswai_aero_s2, &
     6169                     topswad0_aero_s2, solswad0_aero_s2, &
     6170                     topsw_aero_s2, topsw0_aero_s2, &
     6171                     solsw_aero_s2, solsw0_aero_s2, &
     6172                     topswcf_aero_s2, solswcf_aero_s2, &
     6173                                !-C. Kleinschmitt for LW diagnostics
     6174                     toplwad_aero_s2, sollwad_aero_s2,&
     6175                     toplwai_aero_s2, sollwai_aero_s2, &
     6176                     toplwad0_aero_s2, sollwad0_aero_s2,&
     6177                                !-end
     6178                     ZLWFT0_i, ZFLDN0, ZFLUP0, &
     6179                     ZSWFT0_i, ZFSDN0, ZFSUP0, &
     6180                     cloud_cover_sw_s2)
     6181          ENDIF ! ok_3Deffect
     6182#endif
     6183
    59196184       ENDIF ! aerosol_couple
    59206185       itaprad = 0
    … …  
    61406405       d_t_hin(:, :)=0.
    61416406       CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
    6142             dqi0, dqbs0,paprs, 'hin', abortphy,flag_inhib_tend,itap,0 &
     6407            dqi0, dqbs0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0 &
    61436408#ifdef ISO
    61446409     &    ,dxt0,dxtl0,dxti0 &
    … …  
    62636528       
    62646529       SELECT CASE(flag_emit)
    6265        CASE(1) ! emission volc H2O dans LMDZ
     6530       CASE(1) ! emission volc H2O in LMDZ
    62666531          DO ieru=1, nErupt
    62676532             IF (year_cur==year_emit_vol(ieru).AND.&
    … …  
    62716536               
    62726537                IF(flag_verbose_strataer) print *,'IN physiq_mod: date=',year_cur,mth_cur,day_cur
    6273                 ! initialisation tendance q emission
     6538                ! initialisation of q tendency emission
    62746539                d_q_emiss(:,:)=0.
    62756540                ! daily injection mass emission - NL
    … …  
    62786543                !
    62796544                CALL STRATEMIT(pdtphys,pdtphys,latitude_deg,longitude_deg,t_seri,&
    6280                      pplay,paprs,tr_seri,&
    6281                      m_H2O_emiss_vol_daily,&
    6282                      xlat_min_vol(ieru),xlat_max_vol(ieru),&
    6283                      xlon_min_vol(ieru),xlon_max_vol(ieru),&
    6284                      altemiss_vol(ieru),&
    6285                      sigma_alt_vol(ieru),1,&
    6286                      1,nAerErupt+1,0)
     6545                    pplay,paprs,tr_seri,&
     6546                    m_H2O_emiss_vol_daily,&
     6547                    xlat_min_vol(ieru),xlat_max_vol(ieru),&
     6548                    xlon_min_vol(ieru),xlon_max_vol(ieru),&
     6549                    altemiss_vol(ieru),sigma_alt_vol(ieru),1,1.,&
     6550                    nAerErupt+1,0)
    62876551               
    62886552                IF(flag_verbose_strataer) print *,'IN physiq_mod: min max d_q_emiss=',&
    … …  
    62986562    ENDIF
    62996563#endif
    6300 
    63016564
    63026565!===============================================================
    … …  
    67377000    t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
    67387001
    6739     !=======================================================================
    6740     !   SORTIES
    6741     !=======================================================================
    6742     !
    6743     !IM initialisation + calculs divers diag AMIP2
    6744     !
    6745     include "calcul_divers.h"
    6746     !
    6747     !IM Interpolation sur les niveaux de pression du NMC
    6748     !   -------------------------------------------------
    6749     !
    6750     include "calcul_STDlev.h"
    6751     !
    6752     ! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
    6753     CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
    6754     !
     7002    !==================================================================
     7003    !--OB water mass fixer for the physics
     7004    !--water profiles are corrected to force mass conservation of water
     7005    !--currently flag is turned off
     7006    !==================================================================
     7007    IF (mass_fixer) THEN
     7008#ifdef ISO
     7009      CALL abort_gcm('physiq 6936','isos pas prevus dans le mass fixer',1)
     7010      ! Camille Risi mai 2024: on attend d'avoir la 4e dimension qui rendra tout plus simple.
     7011#endif
     7012    qql2(:)=0.0
     7013    DO k = 1, klev
     7014      qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k))*zmasse(:,k)
     7015      IF (nqo >= 3) THEN
     7016        qql2(:)=qql2(:)+qs_seri(:,k)*zmasse(:,k)
     7017      ENDIF
     7018      IF (ok_bs) THEN
     7019        qql2(:)=qql2(:)+qbs_seri(:,k)*zmasse(:,k)
     7020      ENDIF
     7021    ENDDO
     7022
     7023#ifdef CPP_StratAer
     7024    IF (ok_qemiss) THEN
     7025       DO k = 1, klev
     7026          qql1(:) = qql1(:)+d_q_emiss(:,k)*zmasse(:,k)
     7027       ENDDO
     7028    ENDIF
     7029#endif
     7030    IF (ok_qch4) THEN
     7031       DO k = 1, klev
     7032          qql1(:) = qql1(:)+d_q_ch4_dtime(:,k)*zmasse(:,k)
     7033       ENDDO
     7034    ENDIF
     7035   
     7036    DO i = 1, klon
     7037      !--compute ratio of what q+ql should be with conservation to what it is
     7038      IF (ok_bs) THEN
     7039        corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i)-bs_fall(i))*pdtphys)/qql2(i)
     7040      ELSE
     7041        corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
     7042      ENDIF
     7043      DO k = 1, klev
     7044        q_seri(i,k) =q_seri(i,k)*corrqql
     7045        ql_seri(i,k)=ql_seri(i,k)*corrqql
     7046        IF (nqo >= 3) THEN
     7047          qs_seri(i,k)=qs_seri(i,k)*corrqql
     7048        ENDIF
     7049        IF (ok_bs) THEN
     7050          qbs_seri(i,k)=qbs_seri(i,k)*corrqql
     7051        ENDIF
     7052      ENDDO
     7053    ENDDO
     7054    ENDIF
     7055    !--fin mass fixer
     7056
    67557057    !cc prw  = eau precipitable
    67567058    !   prlw = colonne eau liquide
    67577059    !   prlw = colonne eau solide
    67587060    !   prbsw = colonne neige soufflee
     7061    !   water_budget = non-conservation residual from the LMDZ physics
     7062    !                  (should be equal to machine precision if mass fixer is activated)
    67597063    prw(:) = 0.
    67607064    prlw(:) = 0.
    67617065    prsw(:) = 0.
    67627066    prbsw(:) = 0.
     7067    water_budget(:) = 0.0
    67637068    DO k = 1, klev
    67647069       prw(:)  = prw(:)  + q_seri(:,k)*zmasse(:,k)
    67657070       prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
    6766        prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
    6767        prbsw(:)= prbsw(:) + qbs_seri(:,k)*zmasse(:,k)
     7071       water_budget(:) = water_budget(:) + (q_seri(:,k)-qx(:,k,ivap)+ql_seri(:,k)-qx(:,k,iliq))*zmasse(:,k)
     7072       IF (nqo >= 3) THEN
     7073         prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
     7074         water_budget(:) = water_budget(:) + (qs_seri(:,k)-qx(:,k,isol))*zmasse(:,k)
     7075       ENDIF
     7076       IF (nqo >= 4 .AND. ok_bs) THEN
     7077         prbsw(:)= prbsw(:) + qbs_seri(:,k)*zmasse(:,k)
     7078         water_budget(:) = water_budget(:) + (qbs_seri(:,k)-qx(:,k,ibs))*zmasse(:,k)
     7079       ENDIF
    67687080    ENDDO
    6769 
    6770 #ifdef ISO
    6771       DO i = 1, klon
    6772       do ixt=1,ntraciso
    6773        xtprw(ixt,i) = 0.
    6774        DO k = 1, klev
    6775         xtprw(ixt,i) = xtprw(ixt,i) + &
    6776      &           xt_seri(ixt,i,k)*(paprs(i,k)-paprs(i,k+1))/RG
    6777        ENDDO !DO k = 1, klev
    6778       enddo !do ixt=1,ntraciso
    6779       enddo !DO i = 1, klon
    6780 #endif
     7081    water_budget(:)=water_budget(:)+(rain_fall(:)+snow_fall(:)-evap(:))*pdtphys
     7082    IF (ok_bs) THEN
     7083      water_budget(:)=water_budget(:)+bs_fall(:)*pdtphys
     7084    ENDIF
     7085    ! Camille Risi mai 2024: pour les isotopes, on attend d'avoir la 4e dimension, ça rendra tout plus facile
     7086    ! ces variables sont diagnostiques, donc pas indispensables
     7087
     7088    !=======================================================================
     7089    !   SORTIES
     7090    !=======================================================================
     7091    !
     7092    !IM initialisation + calculs divers diag AMIP2
     7093    !
     7094    include "calcul_divers.h"
     7095    !
     7096    !IM Interpolation sur les niveaux de pression du NMC
     7097    !   -------------------------------------------------
     7098    !
     7099    include "calcul_STDlev.h"
     7100    !
     7101    ! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
     7102    CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
     7103    !
    67817104    !
    67827105    IF (ANY(type_trac == ['inca','inco'])) THEN
    … …  
    68817204    !IM global posePB      include "write_bilKP_ave.h"
    68827205    !
    6883 
    6884     !--OB mass fixer
    6885     !--profile is corrected to force mass conservation of water
    6886     IF (mass_fixer) THEN
    6887     qql2(:)=0.0
    6888     DO k = 1, klev
    6889       qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
    6890     ENDDO
    6891    
    6892 #ifdef CPP_StratAer
    6893     IF (ok_qemiss) THEN
    6894        DO k = 1, klev
    6895           qql1(:) = qql1(:)+d_q_emiss(:,k)*zmasse(:,k)
    6896        ENDDO
    6897     ENDIF
    6898 #endif
    6899     IF (ok_qch4) THEN
    6900        DO k = 1, klev
    6901           qql1(:) = qql1(:)+d_q_ch4_dtime(:,k)*zmasse(:,k)
    6902        ENDDO
    6903     ENDIF
    6904    
    6905     DO i = 1, klon
    6906       !--compute ratio of what q+ql should be with conservation to what it is
    6907       corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
    6908       DO k = 1, klev
    6909         q_seri(i,k) =q_seri(i,k)*corrqql
    6910         ql_seri(i,k)=ql_seri(i,k)*corrqql
    6911       ENDDO
    6912     ENDDO
    6913 #ifdef ISO
    6914     do ixt=1,ntraciso
    6915     xtql2(ixt,:)=0.0
    6916     DO k = 1, klev
    6917       xtql2(ixt,:)=xtql2(ixt,:)+(xt_seri(ixt,:,k)+xtl_seri(ixt,:,k)+xts_seri(ixt,:,k))*zmasse(:,k)
    6918     ENDDO
    6919     DO i = 1, klon
    6920       !--compute ratio of what q+ql should be with conservation to what it is
    6921       corrxtql(ixt)=(xtql1(ixt,i)+(xtevap(ixt,i)-xtrain_fall(ixt,i)-xtsnow_fall(ixt,i))*pdtphys)/xtql2(ixt,i)
    6922       DO k = 1, klev
    6923         xt_seri(ixt,i,k) =xt_seri(ixt,i,k)*corrxtql(ixt)
    6924         xtl_seri(ixt,i,k)=xtl_seri(ixt,i,k)*corrxtql(ixt)
    6925       ENDDO
    6926     ENDDO   
    6927     enddo !do ixt=1,ntraciso
    6928 #endif
    6929     ENDIF
    6930     !--fin mass fixer
    6931 
    69327206    ! Sauvegarder les valeurs de t et q a la fin de la physique:
    69337207    !
    … …  
    69447218    xtl_ancien(:,:,:)=xtl_seri(:,:,:)
    69457219    xts_ancien(:,:,:)=xts_seri(:,:,:)
     7220    xtbs_ancien(:,:,:)=xtbs_seri(:,:,:)
    69467221#endif
    69477222    CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
    … …  
    70807355         ok_sync, ptconv, read_climoz, clevSTD,          &
    70817356         ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
    7082          flag_aerosol, flag_aerosol_strat, ok_cdnc,t,u1,v1)
     7357         flag_aerosol, flag_aerosol_strat, ok_cdnc,t, u1, v1)
    70837358#endif
    70847359
    70857360#ifndef CPP_XIOS
    7086     CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
    7087 #endif
    7088 
    7089 #endif
    7090 
    7091 ! Pour XIOS : On remet des variables a .false. apres un premier appel
    7092     IF (debut) THEN
    7093 
    7094       IF (using_xios) THEN
    7095         swaero_diag=.FALSE.
    7096         swaerofree_diag=.FALSE.
    7097         dryaod_diag=.FALSE.
    7098         ok_4xCO2atm= .FALSE.
    7099 !       write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
    7100 
    7101         IF (is_master) THEN
    7102           !--setting up swaero_diag to TRUE in XIOS case
    7103           IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. &
    7104              xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. &
    7105              xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR.  &
    7106                (iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. &
    7107                                    xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0"))))  &
    7108              !!!--for now these fields are not in the XML files so they are omitted
    7109              !!!  xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) &
    7110              swaero_diag=.TRUE.
    7111 
    7112           !--setting up swaerofree_diag to TRUE in XIOS case
    7113           IF (xios_field_is_active("SWdnSFCcleanclr").OR.xios_field_is_active("SWupSFCcleanclr").OR. &
    7114              xios_field_is_active("SWupTOAcleanclr").OR.xios_field_is_active("rsucsaf").OR.   &
    7115              xios_field_is_active("rsdcsaf") .OR. xios_field_is_active("LWdnSFCcleanclr").OR. &
    7116              xios_field_is_active("LWupTOAcleanclr")) &
    7117              swaerofree_diag=.TRUE.
    7118 
    7119           !--setting up dryaod_diag to TRUE in XIOS case
    7120           DO naero = 1, naero_tot-1
    7121            IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
    7122           ENDDO
    7123           !
    7124           !--setting up ok_4xCO2atm to TRUE in XIOS case
    7125           IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. &
    7126              xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
    7127              xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
    7128              xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
    7129              xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
    7130              xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
    7131              ok_4xCO2atm=.TRUE.
    7132         ENDIF
    7133         !$OMP BARRIER
    7134         CALL bcast(swaero_diag)
    7135         CALL bcast(swaerofree_diag)
    7136         CALL bcast(dryaod_diag)
    7137         CALL bcast(ok_4xCO2atm)
    7138 !        write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
    7139       ENDIF !using_xios
    7140     ENDIF
     7361      CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
     7362#endif
     7363
     7364#endif
     7365    ! Petit appelle de sorties pour accompagner le travail sur phyex
     7366    if ( iflag_physiq == 1 ) then
     7367        call output_physiqex(debut,jD_eq,pdtphys,presnivs,paprs,u,v,t,qx,cldfra,0.*t,0.*t,0.*t,pbl_tke,theta)
     7368    endif
    71417369
    71427370    !====================================================================
    … …  
    71827410    ! Disabling calls to the prt_alerte function
    71837411    alert_first_call = .FALSE.
     7412
    71847413   
    71857414    IF (lafin) THEN
    … …  
    71947423         IF (read_climoz >= 1) THEN
    71957424           IF (is_mpi_root) CALL nf95_close(ncid_climoz)
    7196             DEALLOCATE(press_edg_climoz) ! pointer
    7197             DEALLOCATE(press_cen_climoz) ! pointer
     7425            DEALLOCATE(press_edg_climoz)
     7426            DEALLOCATE(press_cen_climoz)
    71987427         ENDIF
    71997428       
    72007429       ENDIF
     7430
     7431       IF (using_xios) THEN
     7432
     7433#ifdef INCA
     7434          IF (type_trac == 'inca') THEN
     7435             IF (is_omp_master .AND. grid_type==unstructured) THEN
     7436                CALL finalize_inca
     7437             ENDIF
     7438          ENDIF
     7439#endif
     7440
     7441          IF (is_omp_master .and. grid_type==unstructured) CALL xios_context_finalize
     7442       ENDIF
     7443
     7444       WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
    72017445       
    7202        IF (using_xios) THEN
    7203          IF (is_omp_master) CALL xios_context_finalize
    7204 
    7205 #ifdef INCA
    7206          if (type_trac == 'inca') then
    7207             IF (is_omp_master .and. grid_type==unstructured) THEN
    7208                CALL finalize_inca
    7209             ENDIF
    7210          endif
    7211 #endif
    7212        ENDIF !using_xios
    7213        WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
    72147446    ENDIF
    72157447

  • LMDZ6/trunk/libf/phylmdiso/reevap.F90

    r4491 r4982  
    1   SUBROUTINE reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
    2    &         d_t_eva,d_q_eva,d_ql_eva,d_qs_eva &
    3 #ifdef ISO
    4              ,xt_seri,xtl_seri,xts_seri,d_xt_eva,d_xtl_eva,d_xts_eva &
    5 #endif
    6    &     )
     1  SUBROUTINE reevap (klon,klev,iflag_ice_thermo,t_seri,qx, &
     2   &         d_t_eva,d_qx_eva)
    73
    84    ! flag to include modifications to ensure energy conservation (if flag >0)
    95    USE add_phys_tend_mod, only : fl_cor_ebil
    106#ifdef ISO
    11     USE infotrac_phy, ONLY: ntiso   
     7    USE infotrac_phy, ONLY: ntiso,nqtot,ivap,iliq,isol,iqWIsoPha
    128#ifdef ISOVERIF
    139    USE isotopes_verif_mod
    … …  
    2319
    2420    INTEGER klon,klev,iflag_ice_thermo
    25     REAL, DIMENSION(klon,klev), INTENT(in) :: t_seri,q_seri,ql_seri,qs_seri
    26     REAL, DIMENSION(klon,klev), INTENT(out) :: d_t_eva,d_q_eva,d_ql_eva,d_qs_eva
     21    REAL, DIMENSION(klon,klev), INTENT(in) :: t_seri
     22    REAL, DIMENSION(klon,klev,nqtot), INTENT(in) ::     qx
     23    REAL, DIMENSION(klon,klev), INTENT(out) :: d_t_eva
     24    REAL, DIMENSION(klon,klev,nqtot), INTENT(out) ::        d_qx_eva
    2725
    2826    REAL za,zb,zdelta,zlvdcp,zlsdcp
    29     INTEGER i,k
     27    INTEGER i,k,ixt,ivapcur,iliqcur,isolcur   
    3028
    31 #ifdef ISO
    32     REAL, DIMENSION(ntiso,klon,klev), INTENT(in) :: xt_seri,xtl_seri,xts_seri
    33     REAL, DIMENSION(ntiso,klon,klev), INTENT(out) :: d_xt_eva,d_xtl_eva,d_xts_eva
    34     integer ixt
    35 #endif
    3629
    3730    !--------Stochastic Boundary Layer Triggering: ALE_BL--------
    … …  
    4235    !IM 100106 BEG : pouvoir sortir les ctes de la physique
    4336    !
     37do ixt=1,1+ntiso
    4438    ! Re-evaporer l'eau liquide nuageuse
    4539    !
     40    iliqcur= iqWIsoPha(ixt,iliq)   
     41    ivapcur= iqWIsoPha(ixt,ivap)   
     42    isolcur= iqWIsoPha(ixt,isol)   
    4643!print *,'rrevap ; fl_cor_ebil:',fl_cor_ebil,' iflag_ice_thermo:',iflag_ice_thermo,' RVTMP2',RVTMP2
    4744    DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
    4845       DO i = 1, klon
    49         if (fl_cor_ebil .GT. 0) then
    50           zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*(q_seri(i,k)+ql_seri(i,k)+qs_seri(i,k)))
    51           zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*(q_seri(i,k)+ql_seri(i,k)+qs_seri(i,k)))
    52         else
    53           zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
    54           !jyg<
    55           !  Attention : Arnaud a propose des formules completement differentes
    56           !                  A verifier !!!
    57           zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
    58         end if
    59           IF (iflag_ice_thermo .EQ. 0) THEN
     46
     47        if (ixt.eq.1) then
     48         if (fl_cor_ebil .GT. 0) then
     49          !zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*(q_seri(i,k)+ql_seri(i,k)+qs_seri(i,k)))
     50          !zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*(q_seri(i,k)+ql_seri(i,k)+qs_seri(i,k)))
     51          zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*(qx(i,k,ivapcur)+qx(i,k,iliqcur)+qx(i,k,isolcur)))
     52          zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*(qx(i,k,ivapcur)+qx(i,k,iliqcur)+qx(i,k,isolcur)))
     53         else
     54           zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*qx(i,k,ivapcur))
     55           !jyg<
     56           !  Attention : Arnaud a propose des formules completement differentes
     57           !                  A verifier !!!
     58           zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*qx(i,k,ivapcur))
     59         end if
     60         IF (iflag_ice_thermo .EQ. 0) THEN
    6061             zlsdcp=zlvdcp
    61           ENDIF
     62         ENDIF
    6263          !>jyg
    63 
    64           IF (iflag_ice_thermo.eq.0) THEN   
    65              !pas necessaire a priori
    66 
    67              zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
    68   zdelta = 0.
    69              zb = MAX(0.0,ql_seri(i,k))
    70              za = - MAX(0.0,ql_seri(i,k)) &
    71                   * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
    72              d_t_eva(i,k) = za
    73              d_q_eva(i,k) = zb
    74              d_ql_eva(i,k) = -ql_seri(i,k)
    75              d_qs_eva(i,k) = 0.
    76 
    77 #ifdef ISO
    78          do ixt=1,ntiso
    79             zb = MAX(0.0,xtl_seri(ixt,i,k))
    80             d_xt_eva(ixt,i,k) = zb
    81             d_xtl_eva(ixt,i,k) = -xtl_seri(ixt,i,k)
    82             d_xts_eva(ixt,i,k) = 0.
    83          enddo
    84 #ifdef ISOVERIF
    85       do ixt=1,ntiso
    86         call iso_verif_noNaN(xt_seri(ixt,i,k), &
    87      &     'reevap 2417: apres evap tot')
    88       enddo
    89       if (iso_eau.gt.0) then
    90               call iso_verif_egalite_choix( &
    91      &           xt_seri(iso_eau,i,k),q_seri(i,k), &
    92      &          'reevap 1891+, après reevap totale',errmax,errmaxrel)
    93               call iso_verif_egalite_choix( &
    94      &           xtl_seri(iso_eau,i,k),ql_seri(i,k), &
    95      &          'reevap 2209+, après reevap totale',errmax,errmaxrel)
    96        endif !if (iso_eau.gt.0) then       
    97       if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then 
    98             if (q_seri(i,k).gt.ridicule) then 
    99                if (iso_verif_o18_aberrant_nostop( &
    100      &              xt_seri(iso_HDO,i,k)/q_seri(i,k), &
    101      &              xt_seri(iso_O18,i,k)/q_seri(i,k), &
    102      &              'reevap 2315: apres reevap totale').eq.1) then
    103                   write(*,*) 'i,k,q_seri(i,k)=',i,k,q_seri(i,k)
    104                   write(*,*) 'd_q_eva(i,k)=',d_q_eva(i,k)
    105                   write(*,*) 'deltaD(d_q_eva(i,k))=',deltaD(d_xt_eva(iso_HDO,i,k)/d_q_eva(i,k))
    106                   write(*,*) 'deltaO18(d_q_eva(i,k))=',deltaO(d_xt_eva(iso_O18,i,k)/d_q_eva(i,k))
    107                   stop
    108               endif !  if (iso_verif_o18_aberrant_nostop
    109             endif !if (q_seri(i,k).gt.errmax) then   
    110         endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then
    111 #ifdef ISOTRAC     
    112              call iso_verif_traceur(xt_seri(1,i,k), &
    113      &           'reevap 2165a')
    114              call iso_verif_traceur_pbidouille(xt_seri(1,i,k), &
    115      &           'reevap 2165b')
    116 #endif               
    117          
    118 #endif           
    119 #endif
    120 
    121           ELSE
    122 
    12364             !CR: on r\'e-\'evapore eau liquide et glace
    12465
    … …  
    12768             !        za = - MAX(0.0,ql_seri(i,k)) &
    12869             !             * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
    129              zb = MAX(0.0,ql_seri(i,k)+qs_seri(i,k))
    130              za = - MAX(0.0,ql_seri(i,k))*zlvdcp &
    131                   - MAX(0.0,qs_seri(i,k))*zlsdcp
     70             za = - MAX(0.0,qx(i,k,iliqcur))*zlvdcp &
     71                  - MAX(0.0,qx(i,k,iliqcur))*zlsdcp
    13272             d_t_eva(i,k) = za
    133              d_q_eva(i,k) = zb
    134              d_ql_eva(i,k) = -ql_seri(i,k)
    135              d_qs_eva(i,k) = -qs_seri(i,k)
    13673
    137 #ifdef ISO
    138          do ixt=1,ntiso
    139             zb = MAX(0.0,xtl_seri(ixt,i,k)+xts_seri(ixt,i,k))
    140             d_xt_eva(ixt,i,k) = zb
    141             d_xtl_eva(ixt,i,k) = -xtl_seri(ixt,i,k)
    142             d_xts_eva(ixt,i,k) = -xts_seri(ixt,i,k)
    143          enddo
     74          endif !if (ixt.eq.1) then
    14475
    145 #ifdef ISOVERIF
    146       do ixt=1,ntiso
    147       call iso_verif_noNaN(xt_seri(ixt,i,k), &
    148      &     'reevap 2417: apres evap tot')
    149       enddo
    150       if (iso_eau.gt.0) then
    151               call iso_verif_egalite_choix( &
    152      &           xt_seri(iso_eau,i,k),q_seri(i,k), &
    153      &          'reevap 1891, après réévap totale',errmax,errmaxrel)
    154               call iso_verif_egalite_choix( &
    155      &           xtl_seri(iso_eau,i,k),ql_seri(i,k), &
    156      &          'reevap 2209, après réévap totale',errmax,errmaxrel)
    157               call iso_verif_egalite_choix( &
    158      &           xts_seri(iso_eau,i,k),qs_seri(i,k), &
    159      &          'reevap 2209b, après réévap totale',errmax,errmaxrel)
    160        endif !if (iso_eau.gt.0) then
    161      
    162       if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then 
    163             if (q_seri(i,k).gt.ridicule) then 
    164                if (iso_verif_o18_aberrant_nostop( &
    165      &              xt_seri(iso_HDO,i,k)/q_seri(i,k), &
    166      &              xt_seri(iso_O18,i,k)/q_seri(i,k), &
    167      &              'reevap 2408: apres reevap totale').eq.1) then
    168                   write(*,*) 'i,k,q_seri(i,k)=',i,k,q_seri(i,k)                       
    169                   stop
    170               endif !  if (iso_verif_o18_aberrant_nostop
    171             endif !if (q_seri(i,k).gt.errmax) then   
    172         endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then       
    173 #ifdef ISOTRAC     
    174              call iso_verif_traceur(xt_seri(1,i,k), &
    175      &           'reevap 2165c')
    176              call iso_verif_traceur_pbidouille(xt_seri(1,i,k), &
    177      &           'reevap 2165d')
    178 #endif                 
    179 #endif                 
    180 #endif
     76            !zb = MAX(0.0,ql_seri(i,k)+qs_seri(i,k))
     77            !d_q_eva(i,k) = zb
     78            !d_ql_eva(i,k) = -ql_seri(i,k)
     79            !d_qs_eva(i,k) = -qs_seri(i,k)
    18180
    182           ENDIF
     81            zb = MAX(0.0,qx(i,k,iliqcur)+qx(i,k,isolcur))
     82            d_qx_eva(i,k,ivapcur) = zb
     83            d_qx_eva(i,k,iliqcur) = -qx(i,k,iliqcur)
     84            d_qx_eva(i,k,isolcur) = -qx(i,k,isolcur)
     85
    18386
    18487       ENDDO
    18588    ENDDO
    18689
     90    enddo ! do ixt=1,1+niso*(nzone +1)
     91   
     92
    18793RETURN
    18894
Note: See TracChangeset for help on using the changeset viewer.