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

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.