source: LMDZ6/trunk/libf/phylmd/physiq_mod.F90 @ 4462

Last change on this file since 4462 was 4458, checked in by evignon, 17 months ago

mise des seuils d'activation des params de SSO sous flag
pour faciliter les tests de sensibilité à venir

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 192.6 KB
RevLine 
[2418]1!
[1279]2! $Id: physiq_mod.F90 4458 2023-03-08 14:54:13Z lguez $
[2418]3!
[1862]4!#define IO_DEBUG
[2418]5MODULE physiq_mod
[766]6
[2469]7  IMPLICIT NONE
[2418]8
9CONTAINS
10
[2469]11  SUBROUTINE physiq (nlon,nlev, &
12       debut,lafin,pdtphys_, &
13       paprs,pplay,pphi,pphis,presnivs, &
14       u,v,rot,t,qx, &
15       flxmass_w, &
16       d_u, d_v, d_t, d_qx, d_ps)
[524]17
[3776]18! For clarity, the "USE" section is now arranged in alphabetical order,
19! with a separate section for CPP keys
20! PLEASE try to follow this rule
21
22    USE ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
23    USE aero_mod
24    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
25  &      fl_ebil, fl_cor_ebil
[3387]26    USE assert_m, only: assert
[3776]27    USE change_srf_frac_mod
28    USE conf_phys_m, only: conf_phys
29    USE carbon_cycle_mod, ONLY : infocfields_init, RCO2_glo, carbon_cycle_rad
30    USE CFMIP_point_locations   ! IM stations CFMIP
31    USE cmp_seri_mod
32    USE dimphy
33    USE etat0_limit_unstruct_mod
34    USE FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
35    USE fonte_neige_mod, ONLY  : fonte_neige_get_vars
[4367]36    USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
[2469]37    USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
38         histwrite, ju2ymds, ymds2ju, getin
[3776]39    USE ioipsl_getin_p_mod, ONLY : getin_p
40    USE indice_sol_mod
[4389]41    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac
[4069]42    USE readTracFiles_mod, ONLY: addPhase
[4056]43    USE strings_mod,  ONLY: strIdx
[3776]44    USE iophy
45    USE limit_read_mod, ONLY : init_limit_read
[3435]46    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, grid1dTo2d_glo, grid_type, unstructured
[3776]47    USE mod_phys_lmdz_mpi_data, only: is_mpi_root
[2469]48    USE mod_phys_lmdz_para
[3776]49    USE netcdf95, only: nf95_close
50    USE netcdf, only: nf90_fill_real     ! IM for NMC files
51    USE open_climoz_m, only: open_climoz ! ozone climatology from a file
52    USE ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
53    USE pbl_surface_mod, ONLY : pbl_surface
54    USE phyaqua_mod, only: zenang_an
[4358]55    USE phyetat0_mod, only: phyetat0
[3776]56    USE phystokenc_mod, ONLY: offline, phystokenc
57    USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, &
[4127]58         year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour, calend
[3776]59!!  USE phys_local_var_mod, ONLY : a long list of variables
60!!              ==> see below, after "CPP Keys" section
61    USE phys_state_var_mod ! Variables sauvegardees de la physique
62    USE phys_output_mod
63    USE phys_output_ctrlout_mod
[3981]64    USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level, &
65         alert_first_call, call_alert, prt_alerte
[3776]66    USE readaerosol_mod, ONLY : init_aero_fromfile
67    USE readaerosolstrato_m, ONLY : init_readaerosolstrato
68    USE radlwsw_m, only: radlwsw
69    USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
70    USE regr_pr_time_av_m, only: regr_pr_time_av
[4367]71    USE surface_data,     ONLY : type_ocean, ok_veget
72    USE time_phylmdz_mod, only: current_time, itau_phy, pdtphys, raz_date, update_time
[3776]73    USE tracinca_mod, ONLY: config_inca
[2788]74    USE tropopause_m,     ONLY: dyn_tropopause
[4059]75    USE ice_sursat_mod,  ONLY: flight_init, airplane
[3776]76    USE vampir
77    USE write_field_phy
[4236]78#ifdef CPP_XIOS
79    USE wxios, ONLY: g_ctx, wxios_set_context
80#endif
[3999]81    USE lscp_mod, ONLY : lscp
[4085]82    USE wake_ini_mod, ONLY : wake_ini
[4448]83    USE yamada_ini_mod, ONLY : yamada_ini
[4449]84    USE atke_turbulence_ini_mod, ONLY : atke_ini
[4089]85    USE thermcell_ini_mod, ONLY : thermcell_ini
[4380]86    USE lscp_ini_mod, ONLY : lscp_ini
[3776]87
88    !USE cmp_seri_mod
89!    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
90!  &      fl_ebil, fl_cor_ebil
91
92!!!!!!!!!!!!!!!!!! "USE" section for CPP keys !!!!!!!!!!!!!!!!!!!!!!!!
93!
94!
[2630]95#ifdef CPP_Dust
[3776]96    USE phytracr_spl_mod, ONLY: phytracr_spl, phytracr_spl_out_init
97    USE phys_output_write_spl_mod
98#else
99    USE phytrac_mod, ONLY : phytrac_init, phytrac
100    USE phys_output_write_mod
[2630]101#endif
[3776]102
103
[4367]104#ifdef INCA
105    USE geometry_mod,      ONLY: longitude, latitude, boundslon, boundslat, ind_cell_glo
106    USE time_phylmdz_mod,  ONLY: ndays
107    USE infotrac_phy,      ONLY: nqCO2
108#endif
[3776]109#ifdef REPROBUS
[4367]110    USE chem_rep, ONLY: Init_chem_rep_xjour, d_q_rep, d_ql_rep, d_qi_rep, &
111                        ptrop, ttrop, ztrop, gravit, itroprep, Z1, Z2, fac, B
[3776]112#endif
[4367]113#if defined INCA || defined REPROBUS
114    USE time_phylmdz_mod,    ONLY: annee_ref, day_ini, day_ref, start_time
115    USE vertical_layers_mod, ONLY: aps, bps, ap, bp
116#endif
[3776]117
118
119#ifdef CPP_RRTM
120    USE YOERAD, ONLY : NRADLP
[4367]121!    USE YOESW, ONLY : RSUN
[3776]122#endif
123
124
[3522]125#ifdef CPP_StratAer
126    USE strataer_mod, ONLY: strataer_init
127#endif
[3776]128
129
130#ifdef CPP_XIOS
[4236]131    USE xios, ONLY: xios_update_calendar, xios_context_finalize
132    USE xios, ONLY: xios_get_field_attr, xios_field_is_active, xios_context
133    USE xios, ONLY: xios_set_current_context
[3776]134    USE wxios, ONLY: missing_val, missing_val_omp
135#endif
136#ifndef CPP_XIOS
137    USE paramLMDZ_phy_mod
138#endif
139!
140!
141!!!!!!!!!!!!!!!!!!  END "USE" for CPP keys !!!!!!!!!!!!!!!!!!!!!!
142
143USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
[2606]144       ! [Variables internes non sauvegardees de la physique]
145       ! Variables locales pour effectuer les appels en serie
[4059]146       t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,tr_seri,rneb_seri, &
[2606]147       ! Dynamic tendencies (diagnostics)
[4059]148       d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn,d_rneb_dyn, &
[2606]149       d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, &
150       ! Physic tendencies
151       d_t_con,d_q_con,d_u_con,d_v_con, &
152       d_tr, &                              !! to be removed?? (jyg)
153       d_t_wake,d_q_wake, &
154       d_t_lwr,d_t_lw0,d_t_swr,d_t_sw0, &
155       d_t_ajsb,d_q_ajsb, &
156       d_t_ajs,d_q_ajs,d_u_ajs,d_v_ajs, &
[4367]157!       d_t_ajs_w,d_q_ajs_w, &
158!       d_t_ajs_x,d_q_ajs_x, &
[2606]159       !
[2705]160       d_t_eva,d_q_eva,d_ql_eva,d_qi_eva, &
[2606]161       d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc, &
162       d_t_lscst,d_q_lscst, &
163       d_t_lscth,d_q_lscth, &
164       plul_st,plul_th, &
165       !
166       d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_t_diss, &
[4056]167       d_t_vdf_x, d_t_vdf_w, &
168       d_q_vdf_x, d_q_vdf_w, &
[2606]169       d_ts, &
170       !
[4367]171!       d_t_oli,d_u_oli,d_v_oli, &
[2606]172       d_t_oro,d_u_oro,d_v_oro, &
[2897]173       d_t_oro_gw,d_u_oro_gw,d_v_oro_gw, &
[2606]174       d_t_lif,d_u_lif,d_v_lif, &
175       d_t_ec, &
176       !
177       du_gwd_hines,dv_gwd_hines,d_t_hin, &
178       dv_gwd_rando,dv_gwd_front, &
179       east_gwstress,west_gwstress, &
180       d_q_ch4, &
181       !  Special RRTM
182       ZLWFT0_i,ZSWFT0_i,ZFLDN0,  &
183       ZFLUP0,ZFSDN0,ZFSUP0,      &
184       !
185       topswad_aero,solswad_aero,   &
186       topswai_aero,solswai_aero,   &
187       topswad0_aero,solswad0_aero, &
188       !LW additional
189       toplwad_aero,sollwad_aero,   &
190       toplwai_aero,sollwai_aero,   &
191       toplwad0_aero,sollwad0_aero, &
192       !
193       topsw_aero,solsw_aero,       &
194       topsw0_aero,solsw0_aero,     &
195       topswcf_aero,solswcf_aero,   &
196       tausum_aero,tau3d_aero,      &
[2854]197       drytausum_aero,              &
[2606]198       !
199       !variables CFMIP2/CMIP5
200       topswad_aerop, solswad_aerop,   &
201       topswai_aerop, solswai_aerop,   &
202       topswad0_aerop, solswad0_aerop, &
203       topsw_aerop, topsw0_aerop,      &
204       solsw_aerop, solsw0_aerop,      &
205       topswcf_aerop, solswcf_aerop,   &
206       !LW diagnostics
207       toplwad_aerop, sollwad_aerop,   &
208       toplwai_aerop, sollwai_aerop,   &
209       toplwad0_aerop, sollwad0_aerop, &
210       !
211       ptstar, pt0, slp, &
212       !
213       bils, &
214       !
215       cldh, cldl,cldm, cldq, cldt,      &
216       JrNt,                             &
217       dthmin, evap, fder, plcl, plfc,   &
218       prw, prlw, prsw,                  &
219       s_lcl, s_pblh, s_pblt, s_therm,   &
220       cdragm, cdragh,                   &
221       zustar, zu10m, zv10m, rh2m, qsat2m, &
[3817]222       zq2m, zt2m, zn2mout, weak_inversion, &
[2606]223       zt2m_min_mon, zt2m_max_mon,   &         ! pour calcul_divers.h
224       t2m_min_mon, t2m_max_mon,  &            ! pour calcul_divers.h
225       !
226       s_pblh_x, s_pblh_w, &
227       s_lcl_x, s_lcl_w,   &
228       !
229       slab_wfbils, tpot, tpote,               &
230       ue, uq, ve, vq, zxffonte,               &
[3257]231       uwat, vwat,                             &
[2606]232       zxfqcalving, zxfluxlat,                 &
233       zxrunofflic,                            &
234       zxtsol, snow_lsc, zxfqfonte, zxqsurf,   &
[3888]235       delta_qsurf,                            &
[2606]236       rain_lsc, rain_num,                     &
237       !
238       sens_x, sens_w, &
239       zxfluxlat_x, zxfluxlat_w, &
240       !
[4009]241       pbl_tke_input, tke_dissip, l_mix, wprime,&
[2606]242       t_therm, q_therm, u_therm, v_therm, &
243       cdragh_x, cdragh_w, &
244       cdragm_x, cdragm_w, &
245       kh, kh_x, kh_w, &
246       !
[2730]247       wake_k, &
[3080]248       alp_wake, &
[2635]249       wake_h, wake_omg, &
250                       ! tendencies of delta T and delta q:
251       d_deltat_wk, d_deltaq_wk, &         ! due to wakes
252       d_deltat_wk_gw, d_deltaq_wk_gw, &   ! due to wake induced gravity waves
253       d_deltat_vdf, d_deltaq_vdf, &       ! due to vertical diffusion
254       d_deltat_the, d_deltaq_the, &       ! due to thermals
255       d_deltat_ajs_cv, d_deltaq_ajs_cv, & ! due to dry adjustment of (w) before convection
256                       ! tendencies of wake fractional area and wake number per unit area:
[3208]257       d_s_wk,  d_dens_a_wk,  d_dens_wk, &  ! due to wakes
258!!!       d_s_vdf, d_dens_a_vdf, d_dens_vdf, & ! due to vertical diffusion
259!!!       d_s_the, d_dens_a_the, d_dens_the, & ! due to thermals
[2635]260       !                                 
[3148]261       ptconv, ratqsc, &
[2824]262       wbeff, convoccur, zmax_th, &
[2606]263       sens, flwp, fiwp,  &
[3080]264       alp_bl_conv,alp_bl_det,  &
[2606]265       alp_bl_fluct_m,alp_bl_fluct_tke,  &
266       alp_bl_stat, n2, s2,  &
267       proba_notrig, random_notrig,  &
[3956]268!!       cv_gen,  &  !moved to phys_state_var_mod
[2606]269       !
[3134]270       dnwd0,  &
271       omega,  &
[2606]272       epmax_diag,  &
[3134]273       !    Deep convective variables used in phytrac
274       pmflxr, pmflxs,  &
[3496]275       wdtrainA, wdtrainS, wdtrainM,  &
[3134]276       upwd, dnwd, &
[2606]277       ep,  &
[3134]278       da, mp, &
279       phi, &
280       wght_cvfd, &
281       phi2, &
282       d1a, dam, &
283       ev, &
284       elij, &
[3496]285       qtaa, &
[3134]286       clw, &
287       epmlmMm, eplaMm, &
288       sij, &
[3387]289       !
[4380]290       rneblsvol, &
291       zqsatl, zqsats, &
292       qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
293       Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
[2606]294       cldemi,  &
295       cldfra, cldtau, fiwc,  &
296       fl, re, flwc,  &
297       ref_liq, ref_ice, theta,  &
298       ref_liq_pi, ref_ice_pi,  &
[3780]299       zphi, zx_rh, zx_rhl, zx_rhi,  &
[2606]300       pmfd, pmfu,  &
301       !
302       t2m, fluxlat,  &
303       fsollw, evap_pot,  &
304       fsolsw, wfbils, wfbilo,  &
[2670]305       wfevap, wfrain, wfsnow,  & 
[3134]306       prfl, psfl, fraca, Vprecip,  &
[2606]307       zw2,  &
[3387]308       !
[2606]309       fluxu, fluxv,  &
310       fluxt,  &
[3387]311       !
[2606]312       uwriteSTD, vwriteSTD, &                !pour calcul_STDlev.h
313       wwriteSTD, phiwriteSTD, &              !pour calcul_STDlev.h
314       qwriteSTD, twriteSTD, rhwriteSTD, &    !pour calcul_STDlev.h
[3387]315       !
[2606]316       beta_prec,  &
317       rneb,  &
[2968]318       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic
[2606]319       !
[2630]320
[3387]321    IMPLICIT NONE
[2469]322    !>======================================================================
323    !!
324    !! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
325    !!
326    !! Objet: Moniteur general de la physique du modele
327    !!AA      Modifications quant aux traceurs :
328    !!AA                  -  uniformisation des parametrisations ds phytrac
329    !!AA                  -  stockage des moyennes des champs necessaires
330    !!AA                     en mode traceur off-line
331    !!======================================================================
332    !!   CLEFS CPP POUR LES IO
333    !!   =====================
[1352]334#define histNMC
[2469]335    !!======================================================================
336    !!    modif   ( P. Le Van ,  12/10/98 )
337    !!
338    !!  Arguments:
339    !!
340    !! nlon----input-I-nombre de points horizontaux
341    !! nlev----input-I-nombre de couches verticales, doit etre egale a klev
342    !! debut---input-L-variable logique indiquant le premier passage
343    !! lafin---input-L-variable logique indiquant le dernier passage
344    !! jD_cur       -R-jour courant a l'appel de la physique (jour julien)
345    !! jH_cur       -R-heure courante a l'appel de la physique (jour julien)
346    !! pdtphys-input-R-pas d'integration pour la physique (seconde)
347    !! paprs---input-R-pression pour chaque inter-couche (en Pa)
348    !! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
349    !! pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
350    !! pphis---input-R-geopotentiel du sol
351    !! presnivs-input_R_pressions approximat. des milieux couches ( en PA)
352    !! u-------input-R-vitesse dans la direction X (de O a E) en m/s
353    !! v-------input-R-vitesse Y (de S a N) en m/s
354    !! t-------input-R-temperature (K)
355    !! qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
356    !! d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
357    !! d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
[2496]358    !! d_ql_dyn-input-R-tendance dynamique pour "ql" (kg/kg/s)
359    !! d_qs_dyn-input-R-tendance dynamique pour "qs" (kg/kg/s)
[2469]360    !! flxmass_w -input-R- flux de masse verticale
361    !! d_u-----output-R-tendance physique de "u" (m/s/s)
362    !! d_v-----output-R-tendance physique de "v" (m/s/s)
363    !! d_t-----output-R-tendance physique de "t" (K/s)
364    !! d_qx----output-R-tendance physique de "qx" (kg/kg/s)
365    !! d_ps----output-R-tendance physique de la pression au sol
366    !!======================================================================
367    integer jjmp1
368    !  parameter (jjmp1=jjm+1-1/jjm) ! => (jjmp1=nbp_lat-1/(nbp_lat-1))
369    !  integer iip1
370    !  parameter (iip1=iim+1)
[782]371
[2469]372    include "regdim.h"
373    include "dimsoil.h"
374    include "clesphys.h"
[4089]375    include "alpale.h"
[3011]376    include "dimpft.h"
[2469]377    !======================================================================
[3479]378    LOGICAL, SAVE :: ok_volcan ! pour activer les diagnostics volcaniques
[3600]379    !$OMP THREADPRIVATE(ok_volcan)
[4056]380    INTEGER, SAVE :: flag_volc_surfstrat ! pour imposer le cool/heat rate à la surf/strato
[3989]381    !$OMP THREADPRIVATE(flag_volc_surfstrat)
[2469]382    LOGICAL ok_cvl  ! pour activer le nouveau driver pour convection KE
383    PARAMETER (ok_cvl=.TRUE.)
384    LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
385    PARAMETER (ok_gust=.FALSE.)
[3600]386    INTEGER, SAVE :: iflag_radia     ! active ou non le rayonnement (MPL)
[2469]387    !$OMP THREADPRIVATE(iflag_radia)
388    !======================================================================
389    LOGICAL check ! Verifier la conservation du modele en eau
390    PARAMETER (check=.FALSE.)
391    LOGICAL ok_stratus ! Ajouter artificiellement les stratus
392    PARAMETER (ok_stratus=.FALSE.)
393    !======================================================================
394    REAL amn, amx
395    INTEGER igout
396    !======================================================================
[3317]397    ! Clef iflag_cycle_diurne controlant l'activation du cycle diurne:
[2469]398    ! en attente du codage des cles par Fred
[3317]399    ! iflag_cycle_diurne est initialise par conf_phys et se trouve
400    ! dans clesphys.h (IM)
[2469]401    !======================================================================
402    ! Modele thermique du sol, a activer pour le cycle diurne:
403    !cc      LOGICAL soil_model
404    !cc      PARAMETER (soil_model=.FALSE.)
405    !======================================================================
406    ! Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
407    ! le calcul du rayonnement est celle apres la precipitation des nuages.
408    ! Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
409    ! la condensation et la precipitation. Cette cle augmente les impacts
410    ! radiatifs des nuages.
411    !cc      LOGICAL new_oliq
412    !cc      PARAMETER (new_oliq=.FALSE.)
413    !======================================================================
414    ! Clefs controlant deux parametrisations de l'orographie:
415    !c      LOGICAL ok_orodr
416    !cc      PARAMETER (ok_orodr=.FALSE.)
417    !cc      LOGICAL ok_orolf
418    !cc      PARAMETER (ok_orolf=.FALSE.)
419    !======================================================================
420    LOGICAL ok_journe ! sortir le fichier journalier
[3461]421    SAVE ok_journe
[2469]422    !$OMP THREADPRIVATE(ok_journe)
423    !
424    LOGICAL ok_mensuel ! sortir le fichier mensuel
[3461]425    SAVE ok_mensuel
[2469]426    !$OMP THREADPRIVATE(ok_mensuel)
427    !
428    LOGICAL ok_instan ! sortir le fichier instantane
[3461]429    SAVE ok_instan
[2469]430    !$OMP THREADPRIVATE(ok_instan)
431    !
432    LOGICAL ok_LES ! sortir le fichier LES
[3461]433    SAVE ok_LES                           
[2469]434    !$OMP THREADPRIVATE(ok_LES)                 
435    !
436    LOGICAL callstats ! sortir le fichier stats
[3461]437    SAVE callstats                           
[2469]438    !$OMP THREADPRIVATE(callstats)                 
439    !
440    LOGICAL ok_region ! sortir le fichier regional
441    PARAMETER (ok_region=.FALSE.)
442    !======================================================================
[3461]443    REAL seuil_inversion
444    SAVE seuil_inversion
[2469]445    !$OMP THREADPRIVATE(seuil_inversion)
[3461]446    INTEGER iflag_ratqs
447    SAVE iflag_ratqs
[2469]448    !$OMP THREADPRIVATE(iflag_ratqs)
449    real facteur
[1507]450
[2469]451    REAL wmax_th(klon)
452    REAL tau_overturning_th(klon)
[878]453
[3461]454    INTEGER lmax_th(klon)
455    INTEGER limbas(klon)
456    REAL ratqscth(klon,klev)
457    REAL ratqsdiff(klon,klev)
458    REAL zqsatth(klon,klev)
[878]459
[2469]460    !======================================================================
461    !
[4143]462    ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional)
463    INTEGER,SAVE :: ivap, iliq, isol, irneb
464!$OMP THREADPRIVATE(ivap, iliq, isol, irneb)
[2469]465    !
466    !
467    ! Variables argument:
468    !
469    INTEGER nlon
470    INTEGER nlev
471    REAL,INTENT(IN) :: pdtphys_
472    ! NB: pdtphys to be used in physics is in time_phylmdz_mod
473    LOGICAL debut, lafin
474    REAL paprs(klon,klev+1)
475    REAL pplay(klon,klev)
476    REAL pphi(klon,klev)
477    REAL pphis(klon)
478    REAL presnivs(klev)
[2799]479!JLD    REAL znivsig(klev)
480!JLD    real pir
[719]481
[2469]482    REAL u(klon,klev)
483    REAL v(klon,klev)
[2333]484
[2469]485    REAL, intent(in):: rot(klon, klev)
486    ! relative vorticity, in s-1, needed for frontal waves
[2333]487
[2469]488    REAL t(klon,klev),thetal(klon,klev)
489    ! thetal: ligne suivante a decommenter si vous avez les fichiers
490    !     MPL 20130625
491    ! fth_fonctions.F90 et parkind1.F90
492    ! sinon thetal=theta
493    !     REAL fth_thetae,fth_thetav,fth_thetal
494    REAL qx(klon,klev,nqtot)
495    REAL flxmass_w(klon,klev)
496    REAL d_u(klon,klev)
497    REAL d_v(klon,klev)
498    REAL d_t(klon,klev)
499    REAL d_qx(klon,klev,nqtot)
500    REAL d_ps(klon)
[2897]501  ! variables pour tend_to_tke
502    REAL duadd(klon,klev)
503    REAL dvadd(klon,klev)
504    REAL dtadd(klon,klev)
505
[2271]506#ifndef CPP_XIOS
[2997]507    REAL, SAVE :: missing_val=nf90_fill_real
[2271]508#endif
[3134]509!!   Variables moved to phys_local_var_mod
510!!    ! Variables pour le transport convectif
511!!    real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
512!!    real wght_cvfd(klon,klev)
513!!    ! Variables pour le lessivage convectif
514!!    ! RomP >>>
515!!    real phi2(klon,klev,klev)
516!!    real d1a(klon,klev),dam(klon,klev)
517!!    real ev(klon,klev)
518!!    real clw(klon,klev),elij(klon,klev,klev)
519!!    real epmlmMm(klon,klev,klev),eplaMm(klon,klev)
520!!    ! RomP <<<
[2469]521    !IM definition dynamique o_trac dans phys_output_open
522    !      type(ctrl_out) :: o_trac(nqtot)
[524]523
[2469]524    ! variables a une pression donnee
525    !
526    include "declare_STDlev.h"
527    !
528    !
529    include "radopt.h"
530    !
531    !
532    INTEGER n
533    !ym      INTEGER npoints
534    !ym      PARAMETER(npoints=klon)
535    !
536    INTEGER nregISCtot
537    PARAMETER(nregISCtot=1)
538    !
539    ! imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties
540    ! sur 1 region rectangulaire y compris pour 1 point
541    ! imin_debut : indice minimum de i; nbpti : nombre de points en
542    ! direction i (longitude)
543    ! jmin_debut : indice minimum de j; nbptj : nombre de points en
544    ! direction j (latitude)
[2799]545!JLD    INTEGER imin_debut, nbpti
546!JLD    INTEGER jmin_debut, nbptj
[2469]547    !IM: region='3d' <==> sorties en global
548    CHARACTER*3 region
549    PARAMETER(region='3d')
[3461]550    LOGICAL ok_hf
[2469]551    !
[3461]552    SAVE ok_hf
[2469]553    !$OMP THREADPRIVATE(ok_hf)
[524]554
[3461]555    INTEGER, PARAMETER :: longcles=20
556    REAL, SAVE :: clesphy0(longcles)
[2469]557    !$OMP THREADPRIVATE(clesphy0)
558    !
559    ! Variables propres a la physique
[3461]560    INTEGER, SAVE :: itap         ! compteur pour la physique
[2469]561    !$OMP THREADPRIVATE(itap)
[2235]562
[2469]563    INTEGER, SAVE :: abortphy=0   ! Reprere si on doit arreter en fin de phys
564    !$OMP THREADPRIVATE(abortphy)
565    !
[3461]566    REAL,SAVE ::  solarlong0
[2469]567    !$OMP THREADPRIVATE(solarlong0)
[987]568
[2469]569    !
570    !  Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
571    !
572    !IM 141004     REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
573    REAL zulow(klon),zvlow(klon)
574    !
575    INTEGER igwd,idx(klon),itest(klon)
576    !
577    !      REAL,allocatable,save :: run_off_lic_0(:)
578    ! !$OMP THREADPRIVATE(run_off_lic_0)
579    !ym      SAVE run_off_lic_0
580    !KE43
581    ! Variables liees a la convection de K. Emanuel (sb):
582    !
[3461]583    REAL, SAVE :: bas, top             ! cloud base and top levels
[2469]584    !$OMP THREADPRIVATE(bas, top)
585    !------------------------------------------------------------------
586    ! Upmost level reached by deep convection and related variable (jyg)
587    !
[4367]588!    INTEGER izero
[2469]589    INTEGER k_upper_cv
590    !------------------------------------------------------------------
[3153]591    ! Compteur de l'occurence de cvpas=1
592    INTEGER Ncvpaseq1
593    SAVE Ncvpaseq1
594    !$OMP THREADPRIVATE(Ncvpaseq1)
[2469]595    !
596    !==========================================================================
597    !CR04.12.07: on ajoute les nouvelles variables du nouveau schema
598    !de convection avec poches froides
599    ! Variables li\'ees \`a la poche froide (jyg)
[879]600
[3496]601!!    REAL mipsh(klon,klev)  ! mass flux shed by the adiab ascent at each level
602!!      Moved to phys_state_var_mod
[2469]603    !
604    REAL wape_prescr, fip_prescr
605    INTEGER it_wape_prescr
606    SAVE wape_prescr, fip_prescr, it_wape_prescr
607    !$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
608    !
609    ! variables supplementaires de concvl
610    REAL Tconv(klon,klev)
[3134]611!!    variable moved to phys_local_var_mod
612!!    REAL sij(klon,klev,klev)
[2812]613!!    !
614!!    ! variables pour tester la conservation de l'energie dans concvl
615!!    REAL, DIMENSION(klon,klev)     :: d_t_con_sat
616!!    REAL, DIMENSION(klon,klev)     :: d_q_con_sat
617!!    REAL, DIMENSION(klon,klev)     :: dql_sat
[970]618
[3461]619    REAL, SAVE :: alp_bl_prescr=0.
620    REAL, SAVE :: ale_bl_prescr=0.
621    REAL, SAVE :: wake_s_min_lsp=0.1
[2469]622    !$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
623    !$OMP THREADPRIVATE(wake_s_min_lsp)
[970]624
[3461]625    REAL ok_wk_lsp(klon)
[1516]626
[2469]627    !RC
628    ! Variables li\'ees \`a la poche froide (jyg et rr)
[879]629
[2635]630    INTEGER,  SAVE               :: iflag_wake_tend  ! wake: if =0, then wake state variables are
631                                                     ! updated within calwake
632    !$OMP THREADPRIVATE(iflag_wake_tend)
[3000]633    INTEGER,  SAVE               :: iflag_alp_wk_cond=0 ! wake: if =0, then Alp_wk is the average lifting
634                                                        ! power provided by the wakes; else, Alp_wk is the
635                                                        ! lifting power conditionned on the presence of a
636                                                        ! gust-front in the grid cell.
637    !$OMP THREADPRIVATE(iflag_alp_wk_cond)
[3489]638
[2635]639    REAL t_w(klon,klev),q_w(klon,klev) ! temperature and moisture profiles in the wake region
640    REAL t_x(klon,klev),q_x(klon,klev) ! temperature and moisture profiles in the off-wake region
[879]641
[2469]642    REAL wake_dth(klon,klev)        ! wake : temp pot difference
[879]643
[2469]644    REAL wake_omgbdth(klon,klev)    ! Wake : flux of Delta_Theta
645    ! transported by LS omega
646    REAL wake_dp_omgb(klon,klev)    ! Wake : vertical gradient of
647    ! large scale omega
648    REAL wake_dtKE(klon,klev)       ! Wake : differential heating
649    ! (wake - unpertubed) CONV
650    REAL wake_dqKE(klon,klev)       ! Wake : differential moistening
651    ! (wake - unpertubed) CONV
652    REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
653    REAL wake_spread(klon,klev)     ! spreading term in wake_delt
654    !
655    !pourquoi y'a pas de save??
656    !
[2730]657!!!    INTEGER, SAVE, DIMENSION(klon)   :: wake_k
658!!!    !$OMP THREADPRIVATE(wake_k)
[2469]659    !
660    !jyg<
661    !cc      REAL wake_pe(klon)              ! Wake potential energy - WAPE
662    !>jyg
[879]663
[3000]664    REAL wake_fip_0(klon)           ! Average Front Incoming Power (unconditionned)
[2469]665    REAL wake_gfl(klon)             ! Gust Front Length
[2635]666!!!    REAL wake_dens(klon)         ! moved to phys_state_var_mod
[2469]667    !
668    !
669    REAL dt_dwn(klon,klev)
670    REAL dq_dwn(klon,klev)
671    REAL M_dwn(klon,klev)
672    REAL M_up(klon,klev)
673    REAL dt_a(klon,klev)
674    REAL dq_a(klon,klev)
675    REAL d_t_adjwk(klon,klev)                !jyg
676    REAL d_q_adjwk(klon,klev)                !jyg
677    LOGICAL,SAVE :: ok_adjwk=.FALSE.
678    !$OMP THREADPRIVATE(ok_adjwk)
[2882]679    INTEGER,SAVE :: iflag_adjwk=0            !jyg
680    !$OMP THREADPRIVATE(iflag_adjwk)         !jyg
[2657]681    REAL,SAVE :: oliqmax=999.,oicemax=999.
682    !$OMP THREADPRIVATE(oliqmax,oicemax)
[2469]683    REAL, SAVE :: alp_offset
684    !$OMP THREADPRIVATE(alp_offset)
[3150]685    REAL, SAVE :: dtcon_multistep_max=1.e6
686    !$OMP THREADPRIVATE(dtcon_multistep_max)
687    REAL, SAVE :: dqcon_multistep_max=1.e6
688    !$OMP THREADPRIVATE(dqcon_multistep_max)
689
[2897]690 
[2469]691    !
692    !RR:fin declarations poches froides
693    !==========================================================================
[1032]694
[2469]695    REAL ztv(klon,klev),ztva(klon,klev)
696    REAL zpspsk(klon,klev)
697    REAL ztla(klon,klev),zqla(klon,klev)
698    REAL zthl(klon,klev)
[1638]699
[2469]700    !cc nrlmd le 10/04/2012
[1638]701
[2469]702    !--------Stochastic Boundary Layer Triggering: ALE_BL--------
703    !---Propri\'et\'es du thermiques au LCL
704    real zlcl_th(klon)          ! Altitude du LCL calcul\'e
705    ! continument (pcon dans
706    ! thermcell_main.F90)
707    real fraca0(klon)           ! Fraction des thermiques au LCL
708    real w0(klon)               ! Vitesse des thermiques au LCL
709    real w_conv(klon)           ! Vitesse verticale de grande \'echelle au LCL
710    real tke0(klon,klev+1)      ! TKE au d\'ebut du pas de temps
711    real therm_tke_max0(klon)   ! TKE dans les thermiques au LCL
712    real env_tke_max0(klon)     ! TKE dans l'environnement au LCL
[1638]713
[2799]714!JLD    !---D\'eclenchement stochastique
715!JLD    integer :: tau_trig(klon)
[1638]716
[2469]717    REAL,SAVE :: random_notrig_max=1.
718    !$OMP THREADPRIVATE(random_notrig_max)
[2294]719
[2469]720    !--------Statistical Boundary Layer Closure: ALP_BL--------
721    !---Profils de TKE dans et hors du thermique
722    real therm_tke_max(klon,klev)   ! Profil de TKE dans les thermiques
723    real env_tke_max(klon,klev)     ! Profil de TKE dans l'environnement
[1638]724
[2897]725    !-------Activer les tendances de TKE due a l'orograp??ie---------
726     INTEGER, SAVE :: addtkeoro
727    !$OMP THREADPRIVATE(addtkeoro)
728     REAL, SAVE :: alphatkeoro
729    !$OMP THREADPRIVATE(alphatkeoro)
730     LOGICAL, SAVE :: smallscales_tkeoro
731    !$OMP THREADPRIVATE(smallscales_tkeoro)
[1638]732
[2897]733
734
[2469]735    !cc fin nrlmd le 10/04/2012
[782]736
[2469]737    ! Variables locales pour la couche limite (al1):
738    !
739    !Al1      REAL pblh(klon)           ! Hauteur de couche limite
740    !Al1      SAVE pblh
741    !34EK
742    !
743    ! Variables locales:
744    !
745    !AA
746    !AA  Pour phytrac
747    REAL u1(klon)             ! vents dans la premiere couche U
748    REAL v1(klon)             ! vents dans la premiere couche V
[524]749
[2469]750    !@$$      LOGICAL offline           ! Controle du stockage ds "physique"
751    !@$$      PARAMETER (offline=.false.)
752    !@$$      INTEGER physid
753    REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
754    REAL frac_nucl(klon,klev) ! idem (nucleation)
755    ! RomP >>>
756    REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
757    ! RomP <<<
[2068]758
[2469]759    !IM cf FH pour Tiedtke 080604
760    REAL rain_tiedtke(klon),snow_tiedtke(klon)
761    !
762    !IM 050204 END
763    REAL devap(klon) ! evaporation et sa derivee
764    REAL dsens(klon) ! chaleur sensible et sa derivee
[1279]765
[2469]766    !
767    ! Conditions aux limites
768    !
769    !
770    REAL :: day_since_equinox
771    ! Date de l'equinoxe de printemps
772    INTEGER, parameter :: mth_eq=3, day_eq=21
773    REAL :: jD_eq
[1279]774
[3461]775    LOGICAL, parameter :: new_orbit = .TRUE.
[524]776
[2469]777    !
778    INTEGER lmt_pas
779    SAVE lmt_pas                ! frequence de mise a jour
780    !$OMP THREADPRIVATE(lmt_pas)
781    real zmasse(klon, nbp_lev),exner(klon, nbp_lev)
782    !     (column-density of mass of air in a cell, in kg m-2)
783    real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
[1797]784
[2469]785    !IM sorties
786    REAL un_jour
787    PARAMETER(un_jour=86400.)
788    INTEGER itapm1 !pas de temps de la physique du(es) mois precedents
789    SAVE itapm1    !mis a jour le dernier pas de temps du mois en cours
790    !$OMP THREADPRIVATE(itapm1)
791    !======================================================================
792    !
793    ! Declaration des procedures appelees
794    !
795    EXTERNAL angle     ! calculer angle zenithal du soleil
796    EXTERNAL alboc     ! calculer l'albedo sur ocean
797    EXTERNAL ajsec     ! ajustement sec
798    EXTERNAL conlmd    ! convection (schema LMD)
799    !KE43
800    EXTERNAL conema3  ! convect4.3
801    EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)
802    !AA
803    ! JBM (3/14) fisrtilp_tr not loaded
804    ! EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
805    !                          ! stockage des coefficients necessaires au
806    !                          ! lessivage OFF-LINE et ON-LINE
807    EXTERNAL hgardfou  ! verifier les temperatures
808    EXTERNAL nuage     ! calculer les proprietes radiatives
809    !C      EXTERNAL o3cm      ! initialiser l'ozone
810    EXTERNAL orbite    ! calculer l'orbite terrestre
811    EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
812    EXTERNAL suphel    ! initialiser certaines constantes
813    EXTERNAL transp    ! transport total de l'eau et de l'energie
814    !IM
815    EXTERNAL haut2bas  !variables de haut en bas
816    EXTERNAL ini_undefSTD  !initialise a 0 une variable a 1 niveau de pression
817    EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
818    !     EXTERNAL moy_undefSTD  !moyenne d'1 var a 1 niveau de pression
819    ! EXTERNAL moyglo_aire
820    ! moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
821    ! par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
822    !
823    !
824    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
825    ! Local variables
826    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
827    !
828    REAL rhcl(klon,klev)    ! humiditi relative ciel clair
829    REAL dialiq(klon,klev)  ! eau liquide nuageuse
830    REAL diafra(klon,klev)  ! fraction nuageuse
[4412]831    REAL radocond(klon,klev)  ! eau condensee nuageuse
[2469]832    !
833    !XXX PB
834    REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
835    !
836    REAL zxfluxt(klon, klev)
837    REAL zxfluxq(klon, klev)
838    REAL zxfluxu(klon, klev)
839    REAL zxfluxv(klon, klev)
[1797]840
[2469]841    ! Le rayonnement n'est pas calcule tous les pas, il faut donc
842    !                      sauvegarder les sorties du rayonnement
843    !ym      SAVE  heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
844    !ym      SAVE  sollwdownclr, toplwdown, toplwdownclr
845    !ym      SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
846    !
847    INTEGER itaprad
848    SAVE itaprad
849    !$OMP THREADPRIVATE(itaprad)
850    !
851    REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
852    REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
853    !
854    REAL zsav_tsol(klon)
855    !
856    REAL dist, rmu0(klon), fract(klon)
857    REAL zrmu0(klon), zfract(klon)
858    REAL zdtime, zdtime1, zdtime2, zlongi
859    !
860    REAL qcheck
861    REAL z_avant(klon), z_apres(klon), z_factor(klon)
862    LOGICAL zx_ajustq
863    !
[2799]864    REAL za
865    REAL zx_t, zx_qs, zdelta, zcor
[2469]866    real zqsat(klon,klev)
867    !
[4367]868    INTEGER i, k, iq, nsrf, l, itr
[2469]869    !
870    REAL t_coup
871    PARAMETER (t_coup=234.0)
[1797]872
[2469]873    !ym A voir plus tard !!
874    !ym      REAL zx_relief(iim,jjmp1)
875    !ym      REAL zx_aire(iim,jjmp1)
876    !
877    ! Grandeurs de sorties
878    REAL s_capCL(klon)
879    REAL s_oliqCL(klon), s_cteiCL(klon)
880    REAL s_trmb1(klon), s_trmb2(klon)
881    REAL s_trmb3(klon)
[2707]882
883    ! La convection n'est pas calculee tous les pas, il faut donc
884    !                      sauvegarder les sorties de la convection
885    !ym      SAVE 
886    !ym      SAVE 
887    !ym      SAVE 
888    !
[2730]889    INTEGER itapcv, itapwk
890    SAVE itapcv, itapwk
891    !$OMP THREADPRIVATE(itapcv, itapwk)
[2707]892
[2469]893    !KE43
894    ! Variables locales pour la convection de K. Emanuel (sb):
[524]895
[2469]896    REAL tvp(klon,klev)       ! virtual temp of lifted parcel
897    CHARACTER*40 capemaxcels  !max(CAPE)
[1412]898
[2469]899    REAL rflag(klon)          ! flag fonctionnement de convect
900    INTEGER iflagctrl(klon)          ! flag fonctionnement de convect
[904]901
[2469]902    ! -- convect43:
903    INTEGER ntra              ! nb traceurs pour convect4.3
904    REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
905    REAL dplcldt(klon), dplcldr(klon)
906    !?     .     condm_con(klon,klev),conda_con(klon,klev),
907    !?     .     mr_con(klon,klev),ep_con(klon,klev)
908    !?     .    ,sadiab(klon,klev),wadiab(klon,klev)
909    ! --
910    !34EK
911    !
912    ! Variables du changement
913    !
914    ! con: convection
915    ! lsc: condensation a grande echelle (Large-Scale-Condensation)
916    ! ajs: ajustement sec
917    ! eva: evaporation de l'eau liquide nuageuse
918    ! vdf: couche limite (Vertical DiFfusion)
[2611]919    !
[2469]920    ! tendance nulles
[2812]921    REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
922    REAL, dimension(klon)     :: dsig0, ddens0
923    INTEGER, dimension(klon)  :: wkoccur1
[2801]924    ! tendance buffer pour appel de add_phys_tend
925    REAL, DIMENSION(klon,klev)  :: d_q_ch4_dtime
[2611]926    !
927    ! Flag pour pouvoir ne pas ajouter les tendances.
928    ! Par defaut, les tendances doivente etre ajoutees et
929    ! flag_inhib_tend = 0
930    ! flag_inhib_tend > 0 : tendances non ajoutees, avec un nombre
931    ! croissant de print quand la valeur du flag augmente
932    !!! attention, ce flag doit etre change avec prudence !!!
933    INTEGER :: flag_inhib_tend = 0 !  0 is the default value
934!!    INTEGER :: flag_inhib_tend = 2
[3134]935    !
936    ! Logical switch to a bug : reseting to 0 convective variables at the
937    ! begining of physiq.
938    LOGICAL, SAVE :: ok_bug_cv_trac = .TRUE.
939    !$OMP THREADPRIVATE(ok_bug_cv_trac)
[3180]940    !
941    ! Logical switch to a bug : changing wake_deltat when thermals are active
942    ! even when there are no wakes.
943    LOGICAL, SAVE :: ok_bug_split_th = .TRUE.
944    !$OMP THREADPRIVATE(ok_bug_split_th)
[524]945
[2469]946    !
947    !********************************************************
948    !     declarations
[524]949
[2469]950    !********************************************************
951    !IM 081204 END
952    !
953    REAL pen_u(klon,klev), pen_d(klon,klev)
954    REAL pde_u(klon,klev), pde_d(klon,klev)
955    INTEGER kcbot(klon), kctop(klon), kdtop(klon)
956    !
[3461]957    REAL ratqsbas,ratqshaut,tau_ratqs
958    SAVE ratqsbas,ratqshaut,tau_ratqs
[2469]959    !$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
[2534]960    REAL, SAVE :: ratqsp0=50000., ratqsdp=20000.
961    !$OMP THREADPRIVATE(ratqsp0, ratqsdp)
[644]962
[2469]963    ! Parametres lies au nouveau schema de nuages (SB, PDF)
[3461]964    REAL, SAVE :: fact_cldcon
965    REAL, SAVE :: facttemps
966    !$OMP THREADPRIVATE(fact_cldcon,facttemps)
967    LOGICAL, SAVE :: ok_newmicro
[2469]968    !$OMP THREADPRIVATE(ok_newmicro)
[524]969
[3461]970    INTEGER, SAVE :: iflag_cld_th
[2469]971    !$OMP THREADPRIVATE(iflag_cld_th)
[2877]972!IM logical ptconv(klon,klev)  !passe dans phys_local_var_mod
[2469]973    !IM cf. AM 081204 BEG
[3461]974    LOGICAL ptconvth(klon,klev)
[3999]975
976    REAL picefra(klon,klev)
[4458]977    REAL zrel_oro(klon)
[2469]978    !IM cf. AM 081204 END
979    !
980    ! Variables liees a l'ecriture de la bande histoire physique
981    !
982    !======================================================================
983    !
984    !
[2799]985!JLD    integer itau_w   ! pas de temps ecriture = itap + itau_phy
[2469]986    !
987    !
988    ! Variables locales pour effectuer les appels en serie
989    !
990    !IM RH a 2m (la surface)
991    REAL Lheat
[524]992
[2469]993    INTEGER        length
994    PARAMETER    ( length = 100 )
995    REAL tabcntr0( length       )
996    !
[2799]997!JLD    INTEGER ndex2d(nbp_lon*nbp_lat)
[2469]998    !IM
999    !
1000    !IM AMIP2 BEG
[2799]1001!JLD    REAL moyglo, mountor
[2469]1002    !IM 141004 BEG
1003    REAL zustrdr(klon), zvstrdr(klon)
1004    REAL zustrli(klon), zvstrli(klon)
1005    REAL zustrph(klon), zvstrph(klon)
1006    REAL aam, torsfc
1007    !IM 141004 END
1008    !IM 190504 BEG
1009    !  INTEGER imp1jmp1
1010    !  PARAMETER(imp1jmp1=(iim+1)*jjmp1)
1011    !ym A voir plus tard
1012    !  REAL zx_tmp((nbp_lon+1)*nbp_lat)
1013    !  REAL airedyn(nbp_lon+1,nbp_lat)
1014    !IM 190504 END
[2799]1015!JLD    LOGICAL ok_msk
1016!JLD    REAL msk(klon)
[2469]1017    !ym A voir plus tard
1018    !ym      REAL zm_wo(jjmp1, klev)
1019    !IM AMIP2 END
1020    !
1021    REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
1022    REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
[2799]1023!JLD    REAL zx_tmp_2d(nbp_lon,nbp_lat)
1024!JLD    REAL zx_lon(nbp_lon,nbp_lat)
1025!JLD    REAL zx_lat(nbp_lon,nbp_lat)
[2469]1026    !
[2630]1027    INTEGER nid_ctesGCM
1028    SAVE nid_ctesGCM
1029    !$OMP THREADPRIVATE(nid_ctesGCM)
[2469]1030    !
1031    !IM 280405 BEG
1032    !  INTEGER nid_bilKPins, nid_bilKPave
1033    !  SAVE nid_bilKPins, nid_bilKPave
1034    !  !$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
1035    !
1036    REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
1037    REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
1038    REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
1039    REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
1040    !
[2799]1041!JLD    REAL zjulian
1042!JLD    SAVE zjulian
1043!JLD!$OMP THREADPRIVATE(zjulian)
[2590]1044
[2799]1045!JLD    INTEGER nhori, nvert
1046!JLD    REAL zsto
1047!JLD    REAL zstophy, zout
[2068]1048
[3981]1049    CHARACTER (LEN=20) :: modname='physiq_mod'
[4056]1050    CHARACTER*80 abort_message
[3461]1051    LOGICAL, SAVE ::  ok_sync, ok_sync_omp
[2469]1052    !$OMP THREADPRIVATE(ok_sync)
[3461]1053    REAL date0
[524]1054
[2469]1055    ! essai writephys
[3461]1056    INTEGER fid_day, fid_mth, fid_ins
1057    PARAMETER (fid_ins = 1, fid_day = 2, fid_mth = 3)
1058    INTEGER prof2d_on, prof3d_on, prof2d_av, prof3d_av
1059    PARAMETER (prof2d_on = 1, prof3d_on = 2, prof2d_av = 3, prof3d_av = 4)
[2469]1060    REAL ztsol(klon)
1061    REAL q2m(klon,nbsrf)  ! humidite a 2m
[524]1062
[2469]1063    !IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
1064    CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
[2799]1065    CHARACTER*40 tinst, tave
[2469]1066    REAL cldtaupi(klon,klev) ! Cloud optical thickness for
1067    ! pre-industrial (pi) aerosols
[524]1068
[2863]1069    INTEGER :: naero
[2469]1070    ! Aerosol optical properties
1071    CHARACTER*4, DIMENSION(naero_grp) :: rfname
1072    REAL, DIMENSION(klon,klev)     :: mass_solu_aero ! total mass
1073    ! concentration
1074    ! for all soluble
1075    ! aerosols[ug/m3]
1076    REAL, DIMENSION(klon,klev)     :: mass_solu_aero_pi
1077    ! - " - (pre-industrial value)
[1279]1078
[2469]1079    ! Parameters
1080    LOGICAL ok_ade, ok_aie    ! Apply aerosol (in)direct effects or not
[2738]1081    LOGICAL ok_alw            ! Apply aerosol LW effect or not
[2469]1082    LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013)
1083    REAL bl95_b0, bl95_b1   ! Parameter in Boucher and Lohmann (1995)
[2738]1084    SAVE ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1
1085    !$OMP THREADPRIVATE(ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1)
[2469]1086    LOGICAL, SAVE :: aerosol_couple ! true  : calcul des aerosols dans INCA
1087    ! false : lecture des aerosol dans un fichier
1088    !$OMP THREADPRIVATE(aerosol_couple)   
[3338]1089    LOGICAL, SAVE :: chemistry_couple ! true  : use INCA chemistry O3
1090    ! false : use offline chemistry O3
1091    !$OMP THREADPRIVATE(chemistry_couple)   
[2469]1092    INTEGER, SAVE :: flag_aerosol
1093    !$OMP THREADPRIVATE(flag_aerosol)
[2644]1094    LOGICAL, SAVE :: flag_bc_internal_mixture
1095    !$OMP THREADPRIVATE(flag_bc_internal_mixture)
[2469]1096    !
1097    !--STRAT AEROSOL
[2530]1098    INTEGER, SAVE :: flag_aerosol_strat
[2469]1099    !$OMP THREADPRIVATE(flag_aerosol_strat)
[3412]1100    !
1101    !--INTERACTIVE AEROSOL FEEDBACK ON RADIATION
1102    LOGICAL, SAVE :: flag_aer_feedback
1103    !$OMP THREADPRIVATE(flag_aer_feedback)
1104
[2469]1105    !c-fin STRAT AEROSOL
1106    !
1107    ! Declaration des constantes et des fonctions thermodynamiques
1108    !
[3461]1109    LOGICAL,SAVE :: first=.TRUE.
[2469]1110    !$OMP THREADPRIVATE(first)
[1279]1111
[2788]1112    ! VARIABLES RELATED TO OZONE CLIMATOLOGIES ; all are OpenMP shared
1113    ! Note that pressure vectors are in Pa and in stricly ascending order
1114    INTEGER,SAVE :: read_climoz                ! Read ozone climatology
[2469]1115    !     (let it keep the default OpenMP shared attribute)
1116    !     Allowed values are 0, 1 and 2
1117    !     0: do not read an ozone climatology
1118    !     1: read a single ozone climatology that will be used day and night
1119    !     2: read two ozone climatologies, the average day and night
1120    !     climatology and the daylight climatology
[2788]1121    INTEGER,SAVE :: ncid_climoz                ! NetCDF file identifier
1122    REAL, POINTER, SAVE :: press_cen_climoz(:) ! Pressure levels
1123    REAL, POINTER, SAVE :: press_edg_climoz(:) ! Edges of pressure intervals
1124    REAL, POINTER, SAVE :: time_climoz(:)      ! Time vector
1125    CHARACTER(LEN=13), PARAMETER :: vars_climoz(2) &
1126                                  = ["tro3         ","tro3_daylight"]
1127    ! vars_climoz(1:read_climoz): variables names in climoz file.
[2819]1128    ! vars_climoz(1:read_climoz-2) if read_climoz>2 (temporary)
1129    REAL :: ro3i ! 0<=ro3i<=360 ; required time index in NetCDF file for
1130                 ! the ozone fields, old method.
[1279]1131
[2469]1132    include "YOMCST.h"
1133    include "YOETHF.h"
1134    include "FCTTRE.h"
1135    !IM 100106 BEG : pouvoir sortir les ctes de la physique
1136    include "conema3.h"
1137    include "fisrtilp.h"
1138    include "nuage.h"
1139    include "compbl.h"
1140    !IM 100106 END : pouvoir sortir les ctes de la physique
1141    !
1142    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1143    ! Declarations pour Simulateur COSP
1144    !============================================================
[3511]1145    real :: mr_ozone(klon,klev), phicosp(klon,klev)
[3370]1146
[2469]1147    !IM stations CFMIP
1148    INTEGER, SAVE :: nCFMIP
1149    !$OMP THREADPRIVATE(nCFMIP)
1150    INTEGER, PARAMETER :: npCFMIP=120
1151    INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
1152    REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
1153    !$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP)
1154    INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
1155    REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
1156    !$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
1157    INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
1158    !$OMP THREADPRIVATE(iGCM, jGCM)
1159    logical, dimension(nfiles)            :: phys_out_filestations
1160    logical, parameter :: lNMC=.FALSE.
[1539]1161
[2469]1162    !IM betaCRF
1163    REAL, SAVE :: pfree, beta_pbl, beta_free
1164    !$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
1165    REAL, SAVE :: lon1_beta,  lon2_beta, lat1_beta, lat2_beta
1166    !$OMP THREADPRIVATE(lon1_beta,  lon2_beta, lat1_beta, lat2_beta)
1167    LOGICAL, SAVE :: mskocean_beta
1168    !$OMP THREADPRIVATE(mskocean_beta)
1169    REAL, dimension(klon, klev) :: beta ! facteur sur cldtaurad et
1170    ! cldemirad pour evaluer les
1171    ! retros liees aux CRF
1172    REAL, dimension(klon, klev) :: cldtaurad   ! epaisseur optique
1173    ! pour radlwsw pour
1174    ! tester "CRF off"
1175    REAL, dimension(klon, klev) :: cldtaupirad ! epaisseur optique
1176    ! pour radlwsw pour
1177    ! tester "CRF off"
1178    REAL, dimension(klon, klev) :: cldemirad   ! emissivite pour
1179    ! radlwsw pour tester
1180    ! "CRF off"
1181    REAL, dimension(klon, klev) :: cldfrarad   ! fraction nuageuse
[1735]1182
[3418]1183#ifdef INCA
[4367]1184    REAL :: calday, zxsnow_dummy(klon)
[3418]1185    ! set de variables utilisees pour l'initialisation des valeurs provenant de INCA
1186    REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_tauinca
1187    REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_pizinca
1188    REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_cginca
1189    REAL, DIMENSION(klon,klev,nbands) :: init_ccminca
1190#endif
1191    REAL, DIMENSION(klon,nbtr) :: init_source
1192
[3048]1193    !lwoff=y : offset LW CRE for radiation code and other schemes
1194    REAL, SAVE :: betalwoff
1195    !OMP THREADPRIVATE(betalwoff)
1196!
[2469]1197    INTEGER :: nbtr_tmp ! Number of tracer inside concvl
1198    REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
[2784]1199    REAL, dimension(klon,klev) :: ch_in ! Condensed humidity entering in phytrac (eau liquide)
[2469]1200    integer iostat
[1539]1201
[4009]1202    REAL, dimension(klon,klev+1) :: tke_dissip_ave, l_mix_ave, wprime_ave
[2469]1203    REAL zzz
1204    !albedo SB >>>
[3461]1205    REAL,DIMENSION(6), SAVE :: SFRWL
[3435]1206!$OMP THREADPRIVATE(SFRWL)
[2469]1207    !albedo SB <<<
[1955]1208
[2485]1209    !--OB variables for mass fixer (hard coded for now)
[3461]1210    LOGICAL, PARAMETER :: mass_fixer=.FALSE.
1211    REAL qql1(klon),qql2(klon),corrqql
[2476]1212
[3110]1213    REAL pi
1214
1215    pi = 4. * ATAN(1.)
1216
[3981]1217    ! set-up call to alerte function
1218    call_alert = (alert_first_call .AND. is_master)
1219   
[2469]1220    ! Ehouarn: set value of jjmp1 since it is no longer a "fixed parameter"
1221    jjmp1=nbp_lat
[2344]1222
[2469]1223    !======================================================================
1224    ! Gestion calendrier : mise a jour du module phys_cal_mod
1225    !
1226    pdtphys=pdtphys_
1227    CALL update_time(pdtphys)
[3435]1228    phys_tstep=NINT(pdtphys)
1229#ifdef CPP_XIOS
[4236]1230! switch to XIOS LMDZ physics context
1231    IF (.NOT. debut .AND. is_omp_master) THEN
1232       CALL wxios_set_context()
1233       CALL xios_update_calendar(itap+1)
1234    ENDIF
[3435]1235#endif
[1355]1236
[2469]1237    !======================================================================
1238    ! Ecriture eventuelle d'un profil verticale en entree de la physique.
1239    ! Utilise notamment en 1D mais peut etre active egalement en 3D
1240    ! en imposant la valeur de igout.
1241    !======================================================================d
[2692]1242    IF (prt_level.ge.1) THEN
[2469]1243       igout=klon/2+1/klon
1244       write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1245       write(lunout,*) 'igout, lat, lon ',igout, latitude_deg(igout), &
1246            longitude_deg(igout)
1247       write(lunout,*) &
1248            'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1249       write(lunout,*) &
1250            nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
[879]1251
[2469]1252       write(lunout,*) 'paprs, play, phi, u, v, t'
[2692]1253       DO k=1,klev
[2469]1254          write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
1255               u(igout,k),v(igout,k),t(igout,k)
[2692]1256       ENDDO
[2469]1257       write(lunout,*) 'ovap (g/kg),  oliq (g/kg)'
[2692]1258       DO k=1,klev
[2469]1259          write(lunout,*) qx(igout,k,1)*1000,qx(igout,k,2)*1000.
[2692]1260       ENDDO
1261    ENDIF
[879]1262
[2769]1263    ! Quick check on pressure levels:
[3461]1264    CALL assert(paprs(:, nbp_lev + 1) < paprs(:, nbp_lev), &
[2769]1265            "physiq_mod paprs bad order")
[879]1266
[2692]1267    IF (first) THEN
[4143]1268       ivap = strIdx(tracers(:)%name, addPhase('H2O', 'g'))
1269       iliq = strIdx(tracers(:)%name, addPhase('H2O', 'l'))
1270       isol = strIdx(tracers(:)%name, addPhase('H2O', 's'))
1271       irneb= strIdx(tracers(:)%name, addPhase('H2O', 'r'))
[3435]1272       CALL init_etat0_limit_unstruct
1273       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
[2469]1274       !CR:nvelles variables convection/poches froides
1275
[3461]1276       WRITE(lunout,*) '================================================='
1277       WRITE(lunout,*) 'Allocation des variables locales et sauvegardees'
1278       WRITE(lunout,*) '================================================='
[2692]1279       CALL phys_local_var_init
[2469]1280       !
1281       !     appel a la lecture du run.def physique
[2692]1282       CALL conf_phys(ok_journe, ok_mensuel, &
[2469]1283            ok_instan, ok_hf, &
1284            ok_LES, &
1285            callstats, &
1286            solarlong0,seuil_inversion, &
1287            fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
1288            iflag_cld_th,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
[3989]1289            ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat, aerosol_couple, &
1290            chemistry_couple, flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
[2644]1291            flag_bc_internal_mixture, bl95_b0, bl95_b1, &
[2469]1292                                ! nv flags pour la convection et les
1293                                ! poches froides
1294            read_climoz, &
1295            alp_offset)
[2692]1296       CALL phys_state_var_init(read_climoz)
1297       CALL phys_output_var_init
[3522]1298       IF (read_climoz>=1 .AND. create_etat0_limit .AND. grid_type==unstructured) &
1299          CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
[3435]1300
[3522]1301#ifdef CPP_StratAer
1302       CALL strataer_init
1303#endif
1304
[2469]1305       print*, '================================================='
1306       !
1307       !CR: check sur le nb de traceurs de l eau
[2692]1308       IF ((iflag_ice_thermo.gt.0).and.(nqo==2)) THEN
[2469]1309          WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
[4120]1310               '(H2O_g, H2O_l, H2O_s) but nqo=', nqo, '. Might as well stop here.'
[3531]1311          abort_message='see above'
1312          CALL abort_physic(modname,abort_message,1)
[2692]1313       ENDIF
[2224]1314
[4062]1315       IF (ok_ice_sursat.AND.(iflag_ice_thermo.EQ.0)) THEN
1316          WRITE (lunout, *) ' ok_ice_sursat=y requires iflag_ice_thermo=1 as well'
[4059]1317          abort_message='see above'
1318          CALL abort_physic(modname,abort_message,1)
1319       ENDIF
1320
[4062]1321       IF (ok_ice_sursat.AND.(nqo.NE.4)) THEN
1322          WRITE (lunout, *) ' ok_ice_sursat=y requires 4 H2O tracers ', &
[4120]1323               '(H2O_g, H2O_l, H2O_s, H2O_r) but nqo=', nqo, '. Might as well stop here.'
[4059]1324          abort_message='see above'
1325          CALL abort_physic(modname,abort_message,1)
1326       ENDIF
1327
[4062]1328       IF (ok_plane_h2o.AND..NOT.ok_ice_sursat) THEN
1329          WRITE (lunout, *) ' ok_plane_h2o=y requires ok_ice_sursat=y '
1330          abort_message='see above'
1331          CALL abort_physic(modname,abort_message,1)
1332       ENDIF
1333
1334       IF (ok_plane_contrail.AND..NOT.ok_ice_sursat) THEN
1335          WRITE (lunout, *) ' ok_plane_contrail=y requires ok_ice_sursat=y '
1336          abort_message='see above'
1337          CALL abort_physic(modname,abort_message,1)
1338       ENDIF
1339
[3154]1340       Ncvpaseq1 = 0
[2469]1341       dnwd0=0.0
1342       ftd=0.0
1343       fqd=0.0
1344       cin=0.
1345       !ym Attention pbase pas initialise dans concvl !!!!
1346       pbase=0
1347       !IM 180608
[904]1348
[2469]1349       itau_con=0
[3461]1350       first=.FALSE.
[1797]1351
[2692]1352    ENDIF  ! first
[1797]1353
[2469]1354    !ym => necessaire pour iflag_con != 2   
1355    pmfd(:,:) = 0.
1356    pen_u(:,:) = 0.
1357    pen_d(:,:) = 0.
1358    pde_d(:,:) = 0.
1359    pde_u(:,:) = 0.
1360    aam=0.
1361    d_t_adjwk(:,:)=0
1362    d_q_adjwk(:,:)=0
[1797]1363
[2469]1364    alp_bl_conv(:)=0.
[2245]1365
[2469]1366    torsfc=0.
1367    forall (k=1: nbp_lev) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
[1797]1368
[644]1369
[2469]1370    IF (debut) THEN
1371       CALL suphel ! initialiser constantes et parametres phys.
[3327]1372! tau_gl : constante de rappel de la temperature a la surface de la glace - en
1373       tau_gl=5.
1374       CALL getin_p('tau_gl', tau_gl)
1375! tau_gl : constante de rappel de la temperature a la surface de la glace - en
1376! secondes
1377       tau_gl=86400.*tau_gl
[3461]1378       WRITE(lunout,*) 'debut physiq_mod tau_gl=',tau_gl
1379
[3000]1380       CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond)
[2469]1381       CALL getin_p('random_notrig_max',random_notrig_max)
[2882]1382       CALL getin_p('ok_adjwk',ok_adjwk)
1383       IF (ok_adjwk) iflag_adjwk=2  ! for compatibility with older versions
1384       ! iflag_adjwk: ! 0 = Default: no convective adjustment of w-region
1385                      ! 1 => convective adjustment but state variables are unchanged
1386                      ! 2 => convective adjustment and state variables are changed
1387       CALL getin_p('iflag_adjwk',iflag_adjwk)
[3150]1388       CALL getin_p('dtcon_multistep_max',dtcon_multistep_max)
1389       CALL getin_p('dqcon_multistep_max',dqcon_multistep_max)
[2613]1390       CALL getin_p('oliqmax',oliqmax)
[2657]1391       CALL getin_p('oicemax',oicemax)
[2534]1392       CALL getin_p('ratqsp0',ratqsp0)
1393       CALL getin_p('ratqsdp',ratqsdp)
[2635]1394       iflag_wake_tend = 0
1395       CALL getin_p('iflag_wake_tend',iflag_wake_tend)
[2799]1396       ok_bad_ecmwf_thermo=.TRUE. ! By default thermodynamical constants are set
1397                                  ! in rrtm/suphec.F90 (and rvtmp2 is set to 0).
1398       CALL getin_p('ok_bad_ecmwf_thermo',ok_bad_ecmwf_thermo)
[3134]1399       CALL getin_p('ok_bug_cv_trac',ok_bug_cv_trac)
[3180]1400       CALL getin_p('ok_bug_split_th',ok_bug_split_th)
[2799]1401       fl_ebil = 0 ! by default, conservation diagnostics are desactivated
1402       CALL getin_p('fl_ebil',fl_ebil)
1403       fl_cor_ebil = 0 ! by default, no correction to ensure energy conservation
1404       CALL getin_p('fl_cor_ebil',fl_cor_ebil)
[2984]1405       iflag_phytrac = 1 ! by default we do want to call phytrac
[2973]1406       CALL getin_p('iflag_phytrac',iflag_phytrac)
[3776]1407#ifdef CPP_Dust
1408       IF (iflag_phytrac.EQ.0) THEN
1409         WRITE(lunout,*) 'In order to run with SPLA, iflag_phytrac will be forced to 1'
1410         iflag_phytrac = 1
1411       ENDIF
[4056]1412#endif
[3011]1413       nvm_lmdz = 13
1414       CALL getin_p('NVM',nvm_lmdz)
[3387]1415
[3461]1416       WRITE(lunout,*) 'iflag_alp_wk_cond=',  iflag_alp_wk_cond
1417       WRITE(lunout,*) 'random_ntrig_max=',   random_notrig_max
1418       WRITE(lunout,*) 'ok_adjwk=',           ok_adjwk
1419       WRITE(lunout,*) 'iflag_adjwk=',        iflag_adjwk
1420       WRITE(lunout,*) 'qtcon_multistep_max=',dtcon_multistep_max
1421       WRITE(lunout,*) 'qdcon_multistep_max=',dqcon_multistep_max
1422       WRITE(lunout,*) 'ratqsp0=',            ratqsp0
1423       WRITE(lunout,*) 'ratqsdp=',            ratqsdp
1424       WRITE(lunout,*) 'iflag_wake_tend=',    iflag_wake_tend
1425       WRITE(lunout,*) 'ok_bad_ecmwf_thermo=',ok_bad_ecmwf_thermo
1426       WRITE(lunout,*) 'ok_bug_cv_trac=',     ok_bug_cv_trac
1427       WRITE(lunout,*) 'ok_bug_split_th=',    ok_bug_split_th
1428       WRITE(lunout,*) 'fl_ebil=',            fl_ebil
1429       WRITE(lunout,*) 'fl_cor_ebil=',        fl_cor_ebil
1430       WRITE(lunout,*) 'iflag_phytrac=',      iflag_phytrac
1431       WRITE(lunout,*) 'NVM=',                nvm_lmdz
1432
[3387]1433       !--PC: defining fields to be exchanged between LMDz, ORCHIDEE and NEMO
1434       WRITE(lunout,*) 'Call to infocfields from physiq'
1435       CALL infocfields_init
1436
[2469]1437    ENDIF
[878]1438
[2692]1439    IF (prt_level.ge.1) print *,'CONVERGENCE PHYSIQUE THERM 1 '
[1279]1440
[2469]1441    !======================================================================
1442    ! Gestion calendrier : mise a jour du module phys_cal_mod
1443    !
1444    !     CALL phys_cal_update(jD_cur,jH_cur)
[1279]1445
[2469]1446    !
1447    ! Si c'est le debut, il faut initialiser plusieurs choses
1448    !          ********
1449    !
1450    IF (debut) THEN
1451       !rv CRinitialisation de wght_th et lalim_conv pour la
1452       !definition de la couche alimentation de la convection a partir
1453       !des caracteristiques du thermique
1454       wght_th(:,:)=1.
1455       lalim_conv(:)=1
1456       !RC
1457       ustar(:,:)=0.
[2569]1458!       u10m(:,:)=0.
1459!       v10m(:,:)=0.
[2469]1460       rain_con(:)=0.
1461       snow_con(:)=0.
1462       topswai(:)=0.
1463       topswad(:)=0.
1464       solswai(:)=0.
1465       solswad(:)=0.
[959]1466
[2469]1467       wmax_th(:)=0.
1468       tau_overturning_th(:)=0.
[645]1469
[4389]1470       IF (ANY(type_trac == ['inca','inco'])) THEN
[2469]1471          ! jg : initialisation jusqu'au ces variables sont dans restart
1472          ccm(:,:,:) = 0.
1473          tau_aero(:,:,:,:) = 0.
1474          piz_aero(:,:,:,:) = 0.
1475          cg_aero(:,:,:,:) = 0.
[2372]1476
[2469]1477          config_inca='none' ! default
1478          CALL getin_p('config_inca',config_inca)
[2372]1479
[2469]1480       ELSE
1481          config_inca='none' ! default
[2692]1482       ENDIF
[782]1483
[3435]1484       tau_aero(:,:,:,:) = 1.e-15
1485       piz_aero(:,:,:,:) = 1.
1486       cg_aero(:,:,:,:)  = 0.
1487
[2469]1488       IF (aerosol_couple .AND. (config_inca /= "aero" &
1489            .AND. config_inca /= "aeNP ")) THEN
1490          abort_message &
1491               = 'if aerosol_couple is activated, config_inca need to be ' &
1492               // 'aero or aeNP'
1493          CALL abort_physic (modname,abort_message,1)
1494       ENDIF
[2443]1495
[2469]1496       rnebcon0(:,:) = 0.0
1497       clwcon0(:,:) = 0.0
1498       rnebcon(:,:) = 0.0
1499       clwcon(:,:) = 0.0
[1863]1500
[2469]1501       !
1502       print*,'iflag_coupl,iflag_clos,iflag_wake', &
1503            iflag_coupl,iflag_clos,iflag_wake
[3317]1504       print*,'iflag_cycle_diurne', iflag_cycle_diurne
[2469]1505       !
1506       IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
1507          abort_message = 'Tiedtke needs iflag_cld_th=-2 or -1'
1508          CALL abort_physic (modname,abort_message,1)
1509       ENDIF
1510       !
1511       !
1512       ! Initialiser les compteurs:
1513       !
1514       itap    = 0
1515       itaprad = 0
[2707]1516       itapcv = 0
[2730]1517       itapwk = 0
[878]1518
[2469]1519       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1520       !! Un petit travail \`a faire ici.
1521       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[878]1522
[2692]1523       IF (iflag_pbl>1) THEN
[2469]1524          PRINT*, "Using method MELLOR&YAMADA"
[2692]1525       ENDIF
[956]1526
[2469]1527       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1528       ! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans
1529       ! phylmd plutot que dyn3d
1530       ! Attention : la version precedente n'etait pas tres propre.
1531       ! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1532       ! pour obtenir le meme resultat.
[2731]1533!jyg for fh<
[3435]1534       WRITE(lunout,*) 'Pas de temps phys_tstep pdtphys ',phys_tstep,pdtphys
1535       IF (abs(phys_tstep-pdtphys)>1.e-10) THEN
[2731]1536          abort_message='pas de temps doit etre entier en seconde pour orchidee et XIOS'
1537          CALL abort_physic(modname,abort_message,1)
1538       ENDIF
1539!>jyg
[3435]1540       IF (MOD(NINT(86400./phys_tstep),nbapp_rad).EQ.0) THEN
1541          radpas = NINT( 86400./phys_tstep)/nbapp_rad
[2469]1542       ELSE
1543          WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1544               'multiple de nbapp_rad'
1545          WRITE(lunout,*) 'changer nbapp_rad ou alors commenter ce test ', &
1546               'mais 1+1<>2'
1547          abort_message='nbre de pas de temps physique n est pas multiple ' &
1548               // 'de nbapp_rad'
[2692]1549          CALL abort_physic(modname,abort_message,1)
[2469]1550       ENDIF
[3435]1551       IF (nbapp_cv .EQ. 0) nbapp_cv=86400./phys_tstep
1552       IF (nbapp_wk .EQ. 0) nbapp_wk=86400./phys_tstep
[2730]1553       print *,'physiq, nbapp_cv, nbapp_wk ',nbapp_cv,nbapp_wk
[3435]1554       IF (MOD(NINT(86400./phys_tstep),nbapp_cv).EQ.0) THEN
1555          cvpas_0 = NINT( 86400./phys_tstep)/nbapp_cv
[3150]1556          cvpas = cvpas_0
[2707]1557       print *,'physiq, cvpas ',cvpas
1558       ELSE
1559          WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1560               'multiple de nbapp_cv'
1561          WRITE(lunout,*) 'changer nbapp_cv ou alors commenter ce test ', &
1562               'mais 1+1<>2'
1563          abort_message='nbre de pas de temps physique n est pas multiple ' &
1564               // 'de nbapp_cv'
[3461]1565          CALL abort_physic(modname,abort_message,1)
[2707]1566       ENDIF
[3435]1567       IF (MOD(NINT(86400./phys_tstep),nbapp_wk).EQ.0) THEN
1568          wkpas = NINT( 86400./phys_tstep)/nbapp_wk
[3457]1569!       print *,'physiq, wkpas ',wkpas
[2730]1570       ELSE
1571          WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1572               'multiple de nbapp_wk'
1573          WRITE(lunout,*) 'changer nbapp_wk ou alors commenter ce test ', &
1574               'mais 1+1<>2'
1575          abort_message='nbre de pas de temps physique n est pas multiple ' &
1576               // 'de nbapp_wk'
[3461]1577          CALL abort_physic(modname,abort_message,1)
[2730]1578       ENDIF
[2469]1579       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[3465]1580       CALL init_iophy_new(latitude_deg,longitude_deg)
[524]1581
[3435]1582          !===================================================================
1583          !IM stations CFMIP
1584          nCFMIP=npCFMIP
1585          OPEN(98,file='npCFMIP_param.data',status='old', &
1586               form='formatted',iostat=iostat)
1587          IF (iostat == 0) THEN
1588             READ(98,*,end=998) nCFMIP
1589998          CONTINUE
1590             CLOSE(98)
1591             CONTINUE
1592             IF(nCFMIP.GT.npCFMIP) THEN
1593                print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
1594                CALL abort_physic("physiq", "", 1)
1595             ELSE
1596                print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
1597             ENDIF
1598
1599             !
1600             ALLOCATE(tabCFMIP(nCFMIP))
1601             ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
1602             ALLOCATE(tabijGCM(nCFMIP))
1603             ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
1604             ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
1605             !
1606             ! lecture des nCFMIP stations CFMIP, de leur numero
1607             ! et des coordonnees geographiques lonCFMIP, latCFMIP
1608             !
1609             CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP,  &
1610                  lonCFMIP, latCFMIP)
1611             !
1612             ! identification des
1613             ! 1) coordonnees lonGCM, latGCM des points CFMIP dans la
1614             ! grille de LMDZ
1615             ! 2) indices points tabijGCM de la grille physique 1d sur
1616             ! klon points
1617             ! 3) indices iGCM, jGCM de la grille physique 2d
1618             !
1619             CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
1620                  tabijGCM, lonGCM, latGCM, iGCM, jGCM)
1621             !
1622          ELSE
1623             ALLOCATE(tabijGCM(0))
1624             ALLOCATE(lonGCM(0), latGCM(0))
1625             ALLOCATE(iGCM(0), jGCM(0))
1626          ENDIF
1627
1628#ifdef CPP_IOIPSL
1629
1630       !$OMP MASTER
1631       ! FH : if ok_sync=.true. , the time axis is written at each time step
1632       ! in the output files. Only at the end in the opposite case
[3461]1633       ok_sync_omp=.FALSE.
[3435]1634       CALL getin('ok_sync',ok_sync_omp)
1635       CALL phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
1636            iGCM,jGCM,lonGCM,latGCM, &
1637            jjmp1,nlevSTD,clevSTD,rlevSTD, phys_tstep,ok_veget, &
1638            type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
1639            ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
1640            read_climoz, phys_out_filestations, &
[3630]1641            aerosol_couple, &
[3435]1642            flag_aerosol_strat, pdtphys, paprs, pphis,  &
1643            pplay, lmax_th, ptconv, ptconvth, ivap,  &
1644            d_u, d_t, qx, d_qx, zmasse, ok_sync_omp)
1645       !$OMP END MASTER
1646       !$OMP BARRIER
1647       ok_sync=ok_sync_omp
1648
1649       freq_outNMC(1) = ecrit_files(7)
1650       freq_outNMC(2) = ecrit_files(8)
1651       freq_outNMC(3) = ecrit_files(9)
1652       WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
1653       WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
1654       WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
1655
1656#ifndef CPP_XIOS
1657       CALL ini_paramLMDZ_phy(phys_tstep,nid_ctesGCM)
1658#endif
1659
1660#endif
1661       ecrit_reg = ecrit_reg * un_jour
1662       ecrit_tra = ecrit_tra * un_jour
1663
1664       !XXXPB Positionner date0 pour initialisation de ORCHIDEE
1665       date0 = jD_ref
1666       WRITE(*,*) 'physiq date0 : ',date0
1667       !
1668
1669!       CALL create_climoz(read_climoz)
[3472]1670      IF (.NOT. create_etat0_limit) CALL init_aero_fromfile(flag_aerosol)  !! initialise aero from file for XIOS interpolation (unstructured_grid)
1671      IF (.NOT. create_etat0_limit) CALL init_readaerosolstrato(flag_aerosol_strat)  !! initialise aero strato from file for XIOS interpolation (unstructured_grid)
[3435]1672
1673#ifdef CPP_COSP
1674      IF (ok_cosp) THEN
[3726]1675!           DO k = 1, klev
1676!             DO i = 1, klon
1677!               phicosp(i,k) = pphi(i,k) + pphis(i)
1678!             ENDDO
1679!           ENDDO
[3465]1680        CALL phys_cosp(itap,phys_tstep,freq_cosp, &
[3435]1681               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1682               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1683               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1684               JrNt,ref_liq,ref_ice, &
1685               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1686               zu10m,zv10m,pphis, &
1687               zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1688               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1689               prfl(:,1:klev),psfl(:,1:klev), &
1690               pmflxr(:,1:klev),pmflxs(:,1:klev), &
1691               mr_ozone,cldtau, cldemi)
[3778]1692      ENDIF
1693#endif
[3512]1694
1695#ifdef CPP_COSP2
1696        IF (ok_cosp) THEN
[3726]1697!           DO k = 1, klev
1698!             DO i = 1, klon
1699!               phicosp(i,k) = pphi(i,k) + pphis(i)
1700!             ENDDO
1701!           ENDDO
[3512]1702          CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
1703               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1704               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1705               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1706               JrNt,ref_liq,ref_ice, &
1707               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1708               zu10m,zv10m,pphis, &
1709               zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1710               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1711               prfl(:,1:klev),psfl(:,1:klev), &
1712               pmflxr(:,1:klev),pmflxs(:,1:klev), &
1713               mr_ozone,cldtau, cldemi)
1714       ENDIF
1715#endif
1716
1717#ifdef CPP_COSPV2
1718        IF (ok_cosp) THEN
1719           DO k = 1, klev
1720             DO i = 1, klon
1721               phicosp(i,k) = pphi(i,k) + pphis(i)
1722             ENDDO
1723           ENDDO
1724          CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
1725               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1726               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1727               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1728               JrNt,ref_liq,ref_ice, &
1729               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1730               zu10m,zv10m,pphis, &
1731               phicosp,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1732               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1733               prfl(:,1:klev),psfl(:,1:klev), &
1734               pmflxr(:,1:klev),pmflxs(:,1:klev), &
1735               mr_ozone,cldtau, cldemi)
1736       ENDIF
1737#endif
1738
[3460]1739       !
[3465]1740       !
1741!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[3460]1742       ! Nouvelle initialisation pour le rayonnement RRTM
[3465]1743!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[3435]1744
[3460]1745       CALL iniradia(klon,klev,paprs(1,1:klev+1))
[4085]1746
1747!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1748       CALL wake_ini(rg,rd,rv,prt_level)
[4448]1749       CALL yamada_ini(klon,lunout,prt_level)
[4449]1750       CALL atke_ini(prt_level, lunout, RG, RD, RPI)
[4089]1751       CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
1752   &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
[4380]1753       IF (ok_new_lscp) then
1754           CALL lscp_ini(pdtphys,ok_ice_sursat)
1755       endif
[4085]1756!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1757
[3956]1758       !
1759!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1760       ! Initialisation des champs dans phytrac* qui sont utilises par phys_output_write*
1761       !
1762!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[3776]1763
1764#ifdef CPP_Dust
1765       ! Quand on utilise SPLA, on force iflag_phytrac=1
1766       CALL phytracr_spl_out_init()
1767       CALL phys_output_write_spl(itap, pdtphys, paprs, pphis,                  &
1768                                pplay, lmax_th, aerosol_couple,                 &
1769                                ok_ade, ok_aie, ivap, ok_sync,                  &
1770                                ptconv, read_climoz, clevSTD,                   &
1771                                ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse,      &
1772                                flag_aerosol, flag_aerosol_strat, ok_cdnc)
1773#else
1774       ! phys_output_write écrit des variables traceurs seulement si iflag_phytrac == 1
1775       ! donc seulement dans ce cas on doit appeler phytrac_init()
[3465]1776       IF (iflag_phytrac == 1 ) THEN
1777          CALL phytrac_init()
[3776]1778       ENDIF
[3465]1779       CALL phys_output_write(itap, pdtphys, paprs, pphis,                    &
1780                              pplay, lmax_th, aerosol_couple,                 &
[3630]1781                              ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ok_sync,&
[3465]1782                              ptconv, read_climoz, clevSTD,                   &
1783                              ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
1784                              flag_aerosol, flag_aerosol_strat, ok_cdnc)
[3776]1785#endif
[3460]1786
[3776]1787
[3435]1788#ifdef CPP_XIOS
1789       IF (is_omp_master) CALL xios_update_calendar(1)
1790#endif
[3465]1791       IF(read_climoz>=1 .AND. create_etat0_limit) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
[3877]1792       CALL create_etat0_limit_unstruct
1793       CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
[3435]1794
[3465]1795!jyg<
[3577]1796       IF (iflag_pbl<=1) THEN
1797          ! No TKE for Standard Physics
1798          pbl_tke(:,:,:)=0.
1799
1800       ELSE IF (klon_glo==1) THEN
[3465]1801          pbl_tke(:,:,is_ave) = 0.
1802          DO nsrf=1,nbsrf
1803            DO k = 1,klev+1
1804                 pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
1805                     +pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
1806            ENDDO
1807          ENDDO
[3988]1808       ELSE
[3465]1809          pbl_tke(:,:,is_ave) = 0. !ym missing init : maybe must be initialized in the same way that for klon_glo==1 ??
1810!>jyg
1811       ENDIF
[2469]1812       !IM begin
1813       print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
1814            ,ratqs(1,1)
1815       !IM end
[878]1816
1817
[2469]1818       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1819       !
1820       ! on remet le calendrier a zero
1821       !
1822       IF (raz_date .eq. 1) THEN
1823          itau_phy = 0
1824       ENDIF
[524]1825
[3435]1826!       IF (ABS(phys_tstep-pdtphys).GT.0.001) THEN
1827!          WRITE(lunout,*) 'Pas physique n est pas correct',phys_tstep, &
1828!               pdtphys
1829!          abort_message='Pas physique n est pas correct '
1830!          !           call abort_physic(modname,abort_message,1)
1831!          phys_tstep=pdtphys
1832!       ENDIF
[2469]1833       IF (nlon .NE. klon) THEN
1834          WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,  &
1835               klon
1836          abort_message='nlon et klon ne sont pas coherents'
[2692]1837          CALL abort_physic(modname,abort_message,1)
[2469]1838       ENDIF
1839       IF (nlev .NE. klev) THEN
1840          WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
1841               klev
1842          abort_message='nlev et klev ne sont pas coherents'
[2692]1843          CALL abort_physic(modname,abort_message,1)
[2469]1844       ENDIF
1845       !
[3435]1846       IF (phys_tstep*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
[2469]1847          WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1848          WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1849          abort_message='Nbre d appels au rayonnement insuffisant'
[2692]1850          CALL abort_physic(modname,abort_message,1)
[2469]1851       ENDIF
[3956]1852
1853!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1854       ! Initialisation pour la convection de K.E. et pour les poches froides
1855       !
1856!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1857
[2469]1858       WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
[3989]1859       WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", ok_cvl
[2469]1860       !
1861       !KE43
1862       ! Initialisation pour la convection de K.E. (sb):
1863       IF (iflag_con.GE.3) THEN
[524]1864
[2469]1865          WRITE(lunout,*)"*** Convection de Kerry Emanuel 4.3  "
1866          WRITE(lunout,*) &
1867               "On va utiliser le melange convectif des traceurs qui"
1868          WRITE(lunout,*)"est calcule dans convect4.3"
1869          WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
[524]1870
[2469]1871          DO i = 1, klon
1872             ema_cbmf(i) = 0.
1873             ema_pcb(i)  = 0.
1874             ema_pct(i)  = 0.
1875             !          ema_workcbmf(i) = 0.
1876          ENDDO
1877          !IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1878          DO i = 1, klon
1879             ibas_con(i) = 1
1880             itop_con(i) = 1
1881          ENDDO
1882          !IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1883          !================================================================
1884          !CR:04.12.07: initialisations poches froides
1885          ! Controle de ALE et ALP pour la fermeture convective (jyg)
[2692]1886          IF (iflag_wake>=1) THEN
[2469]1887             CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
1888                  ,alp_bl_prescr, ale_bl_prescr)
1889             ! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1890             !        print*,'apres ini_wake iflag_cld_th=', iflag_cld_th
[2638]1891             !
1892             ! Initialize tendencies of wake state variables (for some flag values
1893             ! they are not computed).
1894             d_deltat_wk(:,:) = 0.
1895             d_deltaq_wk(:,:) = 0.
1896             d_deltat_wk_gw(:,:) = 0.
1897             d_deltaq_wk_gw(:,:) = 0.
1898             d_deltat_vdf(:,:) = 0.
1899             d_deltaq_vdf(:,:) = 0.
1900             d_deltat_the(:,:) = 0.
1901             d_deltaq_the(:,:) = 0.
1902             d_deltat_ajs_cv(:,:) = 0.
1903             d_deltaq_ajs_cv(:,:) = 0.
1904             d_s_wk(:) = 0.
1905             d_dens_wk(:) = 0.
[3956]1906          ENDIF  !  (iflag_wake>=1)
[973]1907
[2469]1908          !        do i = 1,klon
1909          !           Ale_bl(i)=0.
1910          !           Alp_bl(i)=0.
1911          !        enddo
[1638]1912
[3435]1913       !ELSE
1914       !   ALLOCATE(tabijGCM(0))
1915       !   ALLOCATE(lonGCM(0), latGCM(0))
1916       !   ALLOCATE(iGCM(0), jGCM(0))
[3956]1917       ENDIF  !  (iflag_con.GE.3)
1918       !
[2469]1919       DO i=1,klon
1920          rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1921       ENDDO
[1863]1922
[2469]1923       !34EK
1924       IF (ok_orodr) THEN
[524]1925
[2469]1926          ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1927          ! FH sans doute a enlever de finitivement ou, si on le
1928          ! garde, l'activer justement quand ok_orodr = false.
1929          ! ce rugoro est utilise par la couche limite et fait double emploi
1930          ! avec les param\'etrisations sp\'ecifiques de Francois Lott.
1931          !           DO i=1,klon
1932          !             rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1933          !           ENDDO
1934          ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1935          IF (ok_strato) THEN
1936             CALL SUGWD_strato(klon,klev,paprs,pplay)
1937          ELSE
1938             CALL SUGWD(klon,klev,paprs,pplay)
1939          ENDIF
[1863]1940
[2469]1941          DO i=1,klon
1942             zuthe(i)=0.
1943             zvthe(i)=0.
[2692]1944             IF (zstd(i).gt.10.) THEN
[2469]1945                zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1946                zvthe(i)=(1.-zgam(i))*sin(zthe(i))
[2692]1947             ENDIF
[2469]1948          ENDDO
1949       ENDIF
1950       !
1951       !
[3435]1952       lmt_pas = NINT(86400./phys_tstep * 1.0)   ! tous les jours
[2469]1953       WRITE(lunout,*)'La frequence de lecture surface est de ',  &
1954            lmt_pas
1955       !
1956       capemaxcels = 't_max(X)'
1957       t2mincels = 't_min(X)'
1958       t2maxcels = 't_max(X)'
1959       tinst = 'inst(X)'
1960       tave = 'ave(X)'
1961       !IM cf. AM 081204 BEG
[3465]1962       write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
[2469]1963       !IM cf. AM 081204 END
1964       !
1965       !=============================================================
1966       !   Initialisation des sorties
1967       !=============================================================
1968
[3435]1969#ifdef CPP_XIOS
1970       ! Get "missing_val" value from XML files (from temperature variable)
1971       !$OMP MASTER
1972       CALL xios_get_field_attr("temp",default_value=missing_val_omp)
1973       !$OMP END MASTER
1974       !$OMP BARRIER
1975       missing_val=missing_val_omp
1976#endif
1977
[2679]1978#ifdef CPP_XIOS
[3778]1979! Need to put this initialisation after phyetat0 as in the coupled model the XIOS context is only
1980! initialised at that moment
1981       ! Get "missing_val" value from XML files (from temperature variable)
1982       !$OMP MASTER
1983       CALL xios_get_field_attr("temp",default_value=missing_val_omp)
1984       !$OMP END MASTER
1985       !$OMP BARRIER
1986       missing_val=missing_val_omp
[3988]1987       !
1988       ! Now we activate some double radiation call flags only if some
1989       ! diagnostics are requested, otherwise there is no point in doing this
1990       IF (is_master) THEN
1991         !--setting up swaero_diag to TRUE in XIOS case
1992         IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. &
1993            xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. &
1994            xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR.  &
1995              (iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. &
1996                                  xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0"))))  &
1997            !!!--for now these fields are not in the XML files so they are omitted
1998            !!!  xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) &
1999            swaero_diag=.TRUE.
2000 
2001         !--setting up swaerofree_diag to TRUE in XIOS case
2002         IF (xios_field_is_active("SWdnSFCcleanclr").OR.xios_field_is_active("SWupSFCcleanclr").OR. &
2003            xios_field_is_active("SWupTOAcleanclr").OR.xios_field_is_active("rsucsaf").OR.   &
2004            xios_field_is_active("rsdcsaf") .OR. xios_field_is_active("LWdnSFCcleanclr").OR. &
2005            xios_field_is_active("LWupTOAcleanclr")) &
2006            swaerofree_diag=.TRUE.
2007 
2008         !--setting up dryaod_diag to TRUE in XIOS case
2009         DO naero = 1, naero_tot-1
2010          IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
2011         ENDDO
2012         !
2013         !--setting up ok_4xCO2atm to TRUE in XIOS case
2014         IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. &
2015            xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
2016            xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
2017            xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
2018            xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
2019            xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
2020            ok_4xCO2atm=.TRUE.
2021       ENDIF
2022       !$OMP BARRIER
2023       CALL bcast(swaero_diag)
2024       CALL bcast(swaerofree_diag)
2025       CALL bcast(dryaod_diag)
2026       CALL bcast(ok_4xCO2atm)
[2679]2027#endif
[3988]2028       !
[3435]2029       CALL printflag( tabcntr0,radpas,ok_journe, &
2030            ok_instan, ok_region )
[2469]2031       !
2032       !
2033       ! Prescrire l'ozone dans l'atmosphere
2034       !
2035       !c         DO i = 1, klon
2036       !c         DO k = 1, klev
2037       !c            CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
2038       !c         ENDDO
2039       !c         ENDDO
2040       !
[4389]2041       IF (ANY(type_trac == ['inca','inco'])) THEN ! ModThL
[524]2042#ifdef INCA
[2469]2043          CALL VTe(VTphysiq)
2044          CALL VTb(VTinca)
2045          calday = REAL(days_elapsed) + jH_cur
2046          WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
[959]2047
[4127]2048          call init_const_lmdz( &
2049          ndays, nbsrf, is_oce,is_sic, is_ter,is_lic, calend, &
2050          config_inca)
2051
2052          CALL init_inca_geometry( &
2053               longitude, latitude, &
2054               boundslon, boundslat, &
[4224]2055               cell_area, ind_cell_glo)
[4127]2056
[4224]2057          if (grid_type==unstructured) THEN
2058             CALL chemini(  pplay, &
2059                  nbp_lon, nbp_lat, &
2060                  latitude_deg, &
2061                  longitude_deg, &
2062                  presnivs, &
2063                  calday, &
2064                  klon, &
2065                  nqtot, &
2066                  nqo+nqCO2, &
2067                  pdtphys, &
2068                  annee_ref, &
2069                  year_cur, &
2070                  day_ref,  &
2071                  day_ini, &
2072                  start_time, &
2073                  itau_phy, &
2074                  date0, &
2075                  chemistry_couple, &
2076                  init_source, &
2077                  init_tauinca, &
2078                  init_pizinca, &
2079                  init_cginca, &
2080                  init_ccminca)
2081          ELSE
2082             CALL chemini(  pplay, &
2083                  nbp_lon, nbp_lat, &
2084                  latitude_deg, &
2085                  longitude_deg, &
2086                  presnivs, &
2087                  calday, &
2088                  klon, &
2089                  nqtot, &
2090                  nqo+nqCO2, &
2091                  pdtphys, &
2092                  annee_ref, &
2093                  year_cur, &
2094                  day_ref,  &
2095                  day_ini, &
2096                  start_time, &
2097                  itau_phy, &
2098                  date0, &
2099                  chemistry_couple, &
2100                  init_source, &
2101                  init_tauinca, &
2102                  init_pizinca, &
2103                  init_cginca, &
2104                  init_ccminca, &
2105                  io_lon, &
2106                  io_lat)
2107          ENDIF
[4127]2108
[959]2109
[3418]2110          ! initialisation des variables depuis le restart de inca
2111          ccm(:,:,:) = init_ccminca
2112          tau_aero(:,:,:,:) = init_tauinca
2113          piz_aero(:,:,:,:) = init_pizinca
2114          cg_aero(:,:,:,:) = init_cginca
2115!         
2116
2117
[2469]2118          CALL VTe(VTinca)
2119          CALL VTb(VTphysiq)
[524]2120#endif
[2692]2121       ENDIF
[3988]2122       !
[4389]2123       IF (type_trac == 'repr') THEN
[3666]2124#ifdef REPROBUS
2125          CALL chemini_rep(  &
2126               presnivs, &
2127               pdtphys, &
2128               annee_ref, &
2129               day_ref,  &
2130               day_ini, &
2131               start_time, &
2132               itau_phy, &
2133               io_lon, &
2134               io_lat)
2135#endif
2136       ENDIF
[3465]2137
[2469]2138       !$omp single
[2788]2139       IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz,   &
[2820]2140           press_edg_climoz, time_climoz, ok_daily_climoz, adjust_tropopause)
[2469]2141       !$omp end single
2142       !
2143       !IM betaCRF
2144       pfree=70000. !Pa
2145       beta_pbl=1.
2146       beta_free=1.
2147       lon1_beta=-180.
2148       lon2_beta=+180.
2149       lat1_beta=90.
2150       lat2_beta=-90.
2151       mskocean_beta=.FALSE.
[1539]2152
[2469]2153       !albedo SB >>>
[3461]2154       SELECT CASE(nsw)
2155       CASE(2)
[2469]2156          SFRWL(1)=0.45538747
2157          SFRWL(2)=0.54461211
[3461]2158       CASE(4)
[2469]2159          SFRWL(1)=0.45538747
2160          SFRWL(2)=0.32870591
2161          SFRWL(3)=0.18568763
2162          SFRWL(4)=3.02191470E-02
[3461]2163       CASE(6)
[2469]2164          SFRWL(1)=1.28432794E-03
2165          SFRWL(2)=0.12304168
2166          SFRWL(3)=0.33106142
2167          SFRWL(4)=0.32870591
2168          SFRWL(5)=0.18568763
2169          SFRWL(6)=3.02191470E-02
[3461]2170       END SELECT
[2469]2171       !albedo SB <<<
[2227]2172
[2469]2173       OPEN(99,file='beta_crf.data',status='old', &
2174            form='formatted',err=9999)
2175       READ(99,*,end=9998) pfree
2176       READ(99,*,end=9998) beta_pbl
2177       READ(99,*,end=9998) beta_free
2178       READ(99,*,end=9998) lon1_beta
2179       READ(99,*,end=9998) lon2_beta
2180       READ(99,*,end=9998) lat1_beta
2181       READ(99,*,end=9998) lat2_beta
2182       READ(99,*,end=9998) mskocean_beta
21839998   Continue
2184       CLOSE(99)
21859999   Continue
2186       WRITE(*,*)'pfree=',pfree
2187       WRITE(*,*)'beta_pbl=',beta_pbl
2188       WRITE(*,*)'beta_free=',beta_free
2189       WRITE(*,*)'lon1_beta=',lon1_beta
2190       WRITE(*,*)'lon2_beta=',lon2_beta
2191       WRITE(*,*)'lat1_beta=',lat1_beta
2192       WRITE(*,*)'lat2_beta=',lat2_beta
2193       WRITE(*,*)'mskocean_beta=',mskocean_beta
[3048]2194
2195      !lwoff=y : offset LW CRE for radiation code and other schemes
2196      !lwoff=y : betalwoff=1.
2197      betalwoff=0.
2198      IF (ok_lwoff) THEN
2199         betalwoff=1.
2200      ENDIF
2201      WRITE(*,*)'ok_lwoff=',ok_lwoff
2202      !
2203      !lwoff=y to begin only sollw and sollwdown are set up to CS values
2204      sollw = sollw + betalwoff * (sollw0 - sollw)
2205      sollwdown(:)= sollwdown(:) + betalwoff *(-1.*ZFLDN0(:,1) - &
2206                    sollwdown(:))
[3597]2207
2208
[4085]2209
[2469]2210    ENDIF
2211    !
2212    !   ****************     Fin  de   IF ( debut  )   ***************
2213    !
2214    !
2215    ! Incrementer le compteur de la physique
2216    !
2217    itap   = itap + 1
[2795]2218    IF (is_master .OR. prt_level > 9) THEN
[2783]2219      IF (prt_level > 5 .or. MOD(itap,5) == 0) THEN
[2795]2220         WRITE(LUNOUT,*)'Entering physics elapsed seconds since start ', current_time
2221         WRITE(LUNOUT,100)year_cur,mth_cur,day_cur,hour/3600.
2222 100     FORMAT('Date = ',i4.4,' / ',i2.2, ' / ',i2.2,' : ',f20.17)
[2783]2223      ENDIF
2224    ENDIF
[2469]2225    !
2226    !
2227    ! Update fraction of the sub-surfaces (pctsrf) and
2228    ! initialize, where a new fraction has appeared, all variables depending
2229    ! on the surface fraction.
2230    !
[3435]2231    CALL change_srf_frac(itap, phys_tstep, days_elapsed+1,  &
[2469]2232         pctsrf, fevap, z0m, z0h, agesno,              &
2233         falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
[996]2234
[2469]2235    ! Update time and other variables in Reprobus
[4389]2236    IF (type_trac == 'repr') THEN
[1565]2237#ifdef REPROBUS
[2469]2238       CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
2239       print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
2240       CALL Rtime(debut)
[1565]2241#endif
[2692]2242    ENDIF
[1565]2243
[2469]2244    ! Tendances bidons pour les processus qui n'affectent pas certaines
2245    ! variables.
2246    du0(:,:)=0.
2247    dv0(:,:)=0.
2248    dt0 = 0.
2249    dq0(:,:)=0.
2250    dql0(:,:)=0.
2251    dqi0(:,:)=0.
[2635]2252    dsig0(:) = 0.
2253    ddens0(:) = 0.
2254    wkoccur1(:)=1
[2469]2255    !
2256    ! Mettre a zero des variables de sortie (pour securite)
2257    !
2258    DO i = 1, klon
2259       d_ps(i) = 0.0
2260    ENDDO
2261    DO k = 1, klev
2262       DO i = 1, klon
2263          d_t(i,k) = 0.0
2264          d_u(i,k) = 0.0
2265          d_v(i,k) = 0.0
2266       ENDDO
2267    ENDDO
2268    DO iq = 1, nqtot
2269       DO k = 1, klev
2270          DO i = 1, klon
2271             d_qx(i,k,iq) = 0.0
2272          ENDDO
2273       ENDDO
2274    ENDDO
2275    beta_prec_fisrt(:,:)=0.
2276    beta_prec(:,:)=0.
[3134]2277    !
2278    !   Output variables from the convective scheme should not be set to 0
2279    !   since convection is not always called at every time step.
2280    IF (ok_bug_cv_trac) THEN
2281      da(:,:)=0.
2282      mp(:,:)=0.
2283      phi(:,:,:)=0.
2284      ! RomP >>>
2285      phi2(:,:,:)=0.
2286      epmlmMm(:,:,:)=0.
2287      eplaMm(:,:)=0.
2288      d1a(:,:)=0.
2289      dam(:,:)=0.
2290      pmflxr(:,:)=0.
2291      pmflxs(:,:)=0.
2292      ! RomP <<<
2293    ENDIF
[2469]2294    !
2295    ! Ne pas affecter les valeurs entrees de u, v, h, et q
2296    !
2297    DO k = 1, klev
2298       DO i = 1, klon
2299          t_seri(i,k)  = t(i,k)
2300          u_seri(i,k)  = u(i,k)
2301          v_seri(i,k)  = v(i,k)
2302          q_seri(i,k)  = qx(i,k,ivap)
2303          ql_seri(i,k) = qx(i,k,iliq)
2304          !CR: ATTENTION, on rajoute la variable glace
[4098]2305          IF (nqo.EQ.2) THEN             !--vapour and liquid only
[2469]2306             qs_seri(i,k) = 0.
[4066]2307             rneb_seri(i,k) = 0.
[4098]2308          ELSE IF (nqo.EQ.3) THEN        !--vapour, liquid and ice
[2469]2309             qs_seri(i,k) = qx(i,k,isol)
[4066]2310             rneb_seri(i,k) = 0.
[4098]2311          ELSE IF (nqo.EQ.4) THEN        !--vapour, liquid, ice and rneb
[4059]2312             qs_seri(i,k) = qx(i,k,isol)
2313             rneb_seri(i,k) = qx(i,k,irneb)
[2692]2314          ENDIF
[2469]2315       ENDDO
2316    ENDDO
[2476]2317    !
2318    !--OB mass fixer
2319    IF (mass_fixer) THEN
2320    !--store initial water burden
2321    qql1(:)=0.0
[2499]2322    DO k = 1, klev
2323      qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
[2476]2324    ENDDO
2325    ENDIF
2326    !--fin mass fixer
2327
[2469]2328    tke0(:,:)=pbl_tke(:,:,is_ave)
[4056]2329    IF (nqtot > nqo) THEN
2330       ! water isotopes are not included in tr_seri
2331       itr = 0
2332       DO iq = 1, nqtot
[4071]2333         IF(.NOT.tracers(iq)%isInPhysics) CYCLE
[4056]2334         itr = itr+1
[2469]2335          DO  k = 1, klev
2336             DO  i = 1, klon
[4056]2337                tr_seri(i,k,itr) = qx(i,k,iq)
[2469]2338             ENDDO
2339          ENDDO
2340       ENDDO
2341    ELSE
[4056]2342! DC: make sure the final "1" index was meant for 1st H2O phase (vapor) !!!
[4120]2343       tr_seri(:,:,strIdx(tracers(:)%name,addPhase('H2O','g'))) = 0.0
[2469]2344    ENDIF
[3599]2345!
2346! Temporary solutions adressing ticket #104 and the non initialisation of tr_ancien
2347! LF
2348    IF (debut) THEN
2349      WRITE(lunout,*)' WARNING: tr_ancien initialised to tr_seri'
[4056]2350       itr = 0
2351       do iq = 1, nqtot
[4071]2352         IF(.NOT.tracers(iq)%isInPhysics) CYCLE
[4056]2353         itr = itr+1
2354         tr_ancien(:,:,itr)=tr_seri(:,:,itr)       
2355       enddo
[3599]2356    ENDIF
[2469]2357    !
2358    DO i = 1, klon
2359       ztsol(i) = 0.
2360    ENDDO
2361    DO nsrf = 1, nbsrf
2362       DO i = 1, klon
2363          ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
2364       ENDDO
2365    ENDDO
[2611]2366    ! Initialize variables used for diagnostic purpose
[2692]2367    IF (flag_inhib_tend .ne. 0) CALL init_cmp_seri
[524]2368
[2469]2369    ! Diagnostiquer la tendance dynamique
2370    !
2371    IF (ancien_ok) THEN
[2499]2372    !
[3435]2373       d_u_dyn(:,:)  = (u_seri(:,:)-u_ancien(:,:))/phys_tstep
2374       d_v_dyn(:,:)  = (v_seri(:,:)-v_ancien(:,:))/phys_tstep
2375       d_t_dyn(:,:)  = (t_seri(:,:)-t_ancien(:,:))/phys_tstep
2376       d_q_dyn(:,:)  = (q_seri(:,:)-q_ancien(:,:))/phys_tstep
2377       d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
2378       d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
[2499]2379       CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
[3435]2380       d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
[2499]2381       CALL water_int(klon,klev,ql_seri,zmasse,zx_tmp_fi2d)
[3435]2382       d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/phys_tstep
[2499]2383       CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
[3435]2384       d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
[2469]2385       ! !! RomP >>>   td dyn traceur
[4056]2386       IF (nqtot > nqo) d_tr_dyn(:,:,:)=(tr_seri(:,:,:)-tr_ancien(:,:,:))/phys_tstep
[2469]2387       ! !! RomP <<<
[4059]2388       !!d_rneb_dyn(:,:)=(rneb_seri(:,:)-rneb_ancien(:,:))/phys_tstep
2389       d_rneb_dyn(:,:)=0.0
[2469]2390    ELSE
[2499]2391       d_u_dyn(:,:)  = 0.0
2392       d_v_dyn(:,:)  = 0.0
2393       d_t_dyn(:,:)  = 0.0
2394       d_q_dyn(:,:)  = 0.0
2395       d_ql_dyn(:,:) = 0.0
2396       d_qs_dyn(:,:) = 0.0
2397       d_q_dyn2d(:)  = 0.0
2398       d_ql_dyn2d(:) = 0.0
2399       d_qs_dyn2d(:) = 0.0
[2469]2400       ! !! RomP >>>   td dyn traceur
[4056]2401       IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
[2469]2402       ! !! RomP <<<
[4059]2403       d_rneb_dyn(:,:)=0.0
[2469]2404       ancien_ok = .TRUE.
2405    ENDIF
2406    !
2407    ! Ajouter le geopotentiel du sol:
2408    !
2409    DO k = 1, klev
2410       DO i = 1, klon
2411          zphi(i,k) = pphi(i,k) + pphis(i)
2412       ENDDO
2413    ENDDO
2414    !
2415    ! Verifier les temperatures
2416    !
2417    !IM BEG
2418    IF (check) THEN
2419       amn=MIN(ftsol(1,is_ter),1000.)
2420       amx=MAX(ftsol(1,is_ter),-1000.)
2421       DO i=2, klon
2422          amn=MIN(ftsol(i,is_ter),amn)
2423          amx=MAX(ftsol(i,is_ter),amx)
2424       ENDDO
2425       !
2426       PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
2427    ENDIF !(check) THEN
2428    !IM END
2429    !
2430    CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
2431    IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
[2235]2432
[2469]2433    !
2434    !IM BEG
2435    IF (check) THEN
2436       amn=MIN(ftsol(1,is_ter),1000.)
2437       amx=MAX(ftsol(1,is_ter),-1000.)
2438       DO i=2, klon
2439          amn=MIN(ftsol(i,is_ter),amn)
2440          amx=MAX(ftsol(i,is_ter),amx)
2441       ENDDO
2442       !
2443       PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
2444    ENDIF !(check) THEN
2445    !IM END
2446    !
2447    ! Mettre en action les conditions aux limites (albedo, sst, etc.).
2448    ! Prescrire l'ozone et calculer l'albedo sur l'ocean.
2449    !
[2661]2450    ! Update ozone if day change
2451    IF (MOD(itap-1,lmt_pas) == 0) THEN
[2774]2452       IF (read_climoz <= 0) THEN
2453          ! Once per day, update ozone from Royer:
2454          IF (solarlong0<-999.) then
2455             ! Generic case with evolvoing season
2456             zzz=real(days_elapsed+1)
2457          ELSE IF (abs(solarlong0-1000.)<1.e-4) then
2458             ! Particular case with annual mean insolation
2459             zzz=real(90) ! could be revisited
2460             IF (read_climoz/=-1) THEN
2461                abort_message ='read_climoz=-1 is recommended when ' &
2462                     // 'solarlong0=1000.'
2463                CALL abort_physic (modname,abort_message,1)
2464             ENDIF
2465          ELSE
2466             ! Case where the season is imposed with solarlong0
2467             zzz=real(90) ! could be revisited
2468          ENDIF
[2661]2469
[2774]2470          wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
[3666]2471#ifdef REPROBUS
2472          ptrop=dyn_tropopause(t_seri, ztsol, paprs, pplay, rot)/100.
2473          DO i = 1, klon
2474             Z1=t_seri(i,itroprep(i)+1)
2475             Z2=t_seri(i,itroprep(i))
2476             fac=(Z1-Z2)/alog(pplay(i,itroprep(i)+1)/pplay(i,itroprep(i)))
2477             B=Z2-fac*alog(pplay(i,itroprep(i)))
2478             ttrop(i)= fac*alog(ptrop(i))+B
2479!       
2480             Z1= 1.e-3 * ( pphi(i,itroprep(i)+1)+pphis(i) ) / gravit
2481             Z2= 1.e-3 * ( pphi(i,itroprep(i))  +pphis(i) ) / gravit
2482             fac=(Z1-Z2)/alog(pplay(i,itroprep(i)+1)/pplay(i,itroprep(i)))
2483             B=Z2-fac*alog(pplay(i,itroprep(i)))
2484             ztrop(i)=fac*alog(ptrop(i))+B
2485          ENDDO
2486#endif
[2774]2487       ELSE
[2820]2488          !--- ro3i = elapsed days number since current year 1st january, 0h
2489          ro3i=days_elapsed+jh_cur-jh_1jan
2490          !--- scaling for old style files (360 records)
2491          IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
[2788]2492          IF(adjust_tropopause) THEN
[2820]2493             CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),   &
[3086]2494                      ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1),    &
2495                      time_climoz ,  longitude_deg,   latitude_deg,          &
[2968]2496                      dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
[2774]2497          ELSE
[3086]2498             CALL regr_pr_time_av(ncid_climoz,  vars_climoz(1:read_climoz),  &
2499                      ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1),    &
2500                      time_climoz )
[3461]2501          ENDIF
[2774]2502          ! Convert from mole fraction of ozone to column density of ozone in a
2503          ! cell, in kDU:
2504          FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
2505               * zmasse / dobson_u / 1e3
[2788]2506          ! (By regridding ozone values for LMDZ only once a day, we
[2774]2507          ! have already neglected the variation of pressure in one
2508          ! day. So do not recompute "wo" at each time step even if
2509          ! "zmasse" changes a little.)
2510       ENDIF
[2469]2511    ENDIF
2512    !
2513    ! Re-evaporer l'eau liquide nuageuse
2514    !
[2705]2515     CALL reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
2516   &         d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
[1849]2517
[2705]2518     CALL add_phys_tend &
2519            (du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
[2812]2520               'eva',abortphy,flag_inhib_tend,itap,0)
[3461]2521    CALL prt_enerbil('eva',itap)
[2086]2522
[2469]2523    !=========================================================================
2524    ! Calculs de l'orbite.
2525    ! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2526    ! doit donc etre plac\'e avant radlwsw et pbl_surface
[883]2527
[2469]2528    ! !!   jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2692]2529    CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
[2469]2530    day_since_equinox = (jD_cur + jH_cur) - jD_eq
2531    !
2532    !   choix entre calcul de la longitude solaire vraie ou valeur fixee a
2533    !   solarlong0
[2692]2534    IF (solarlong0<-999.) THEN
2535       IF (new_orbit) THEN
[2469]2536          ! calcul selon la routine utilisee pour les planetes
[2692]2537          CALL solarlong(day_since_equinox, zlongi, dist)
2538       ELSE
[2469]2539          ! calcul selon la routine utilisee pour l'AR4
2540          CALL orbite(REAL(days_elapsed+1),zlongi,dist)
[2692]2541       ENDIF
2542    ELSE
[2469]2543       zlongi=solarlong0  ! longitude solaire vraie
2544       dist=1.            ! distance au soleil / moyenne
[2692]2545    ENDIF
[1529]2546
[2692]2547    IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
[524]2548
[2692]2549
[2469]2550    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2551    ! Calcul de l'ensoleillement :
2552    ! ============================
2553    ! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2554    ! l'annee a partir d'une formule analytique.
2555    ! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
2556    ! non nul aux poles.
[2692]2557    IF (abs(solarlong0-1000.)<1.e-4) THEN
2558       CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
[2469]2559            latitude_deg,longitude_deg,rmu0,fract)
[2979]2560       swradcorr(:) = 1.0
2561       JrNt(:) = 1.0
2562       zrmu0(:) = rmu0(:)
[2469]2563    ELSE
2564       ! recode par Olivier Boucher en sept 2015
2565       SELECT CASE (iflag_cycle_diurne)
2566       CASE(0) 
2567          !  Sans cycle diurne
2568          CALL angle(zlongi, latitude_deg, fract, rmu0)
2569          swradcorr = 1.0
2570          JrNt = 1.0
2571          zrmu0 = rmu0
2572       CASE(1) 
2573          !  Avec cycle diurne sans application des poids
2574          !  bit comparable a l ancienne formulation cycle_diurne=true
2575          !  on integre entre gmtime et gmtime+radpas
[3435]2576          zdtime=phys_tstep*REAL(radpas) ! pas de temps du rayonnement (s)
[2469]2577          CALL zenang(zlongi,jH_cur,0.0,zdtime, &
2578               latitude_deg,longitude_deg,rmu0,fract)
2579          zrmu0 = rmu0
2580          swradcorr = 1.0
2581          ! Calcul du flag jour-nuit
2582          JrNt = 0.0
2583          WHERE (fract.GT.0.0) JrNt = 1.0
2584       CASE(2) 
2585          !  Avec cycle diurne sans application des poids
2586          !  On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
2587          !  Comme cette routine est appele a tous les pas de temps de
2588          !  la physique meme si le rayonnement n'est pas appele je
2589          !  remonte en arriere les radpas-1 pas de temps
2590          !  suivant. Petite ruse avec MOD pour prendre en compte le
2591          !  premier pas de temps de la physique pendant lequel
2592          !  itaprad=0
[3435]2593          zdtime1=phys_tstep*REAL(-MOD(itaprad,radpas)-1)     
2594          zdtime2=phys_tstep*REAL(radpas-MOD(itaprad,radpas)-1)
[2469]2595          CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2596               latitude_deg,longitude_deg,rmu0,fract)
2597          !
2598          ! Calcul des poids
2599          !
[3435]2600          zdtime1=-phys_tstep !--on corrige le rayonnement pour representer le
[2469]2601          zdtime2=0.0    !--pas de temps de la physique qui se termine
2602          CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2603               latitude_deg,longitude_deg,zrmu0,zfract)
2604          swradcorr = 0.0
2605          WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &
2606               swradcorr=zfract/fract*zrmu0/rmu0
2607          ! Calcul du flag jour-nuit
2608          JrNt = 0.0
2609          WHERE (zfract.GT.0.0) JrNt = 1.0
2610       END SELECT
2611    ENDIF
[3110]2612    sza_o = ACOS (rmu0) *180./pi
[782]2613
[2692]2614    IF (mydebug) THEN
2615       CALL writefield_phy('u_seri',u_seri,nbp_lev)
2616       CALL writefield_phy('v_seri',v_seri,nbp_lev)
2617       CALL writefield_phy('t_seri',t_seri,nbp_lev)
2618       CALL writefield_phy('q_seri',q_seri,nbp_lev)
2619    ENDIF
[883]2620
[2469]2621    !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2622    ! Appel au pbl_surface : Planetary Boudary Layer et Surface
2623    ! Cela implique tous les interactions des sous-surfaces et la
2624    ! partie diffusion turbulent du couche limit.
2625    !
2626    ! Certains varibales de sorties de pbl_surface sont utiliser que pour
2627    ! ecriture des fihiers hist_XXXX.nc, ces sont :
2628    !   qsol,      zq2m,      s_pblh,  s_lcl,
2629    !   s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
2630    !   s_therm,   s_trmb1,   s_trmb2, s_trmb3,
2631    !   zu10m,     zv10m,   fder,
[3888]2632    !   zxqsurf,   delta_qsurf,
2633    !   rh2m,      zxfluxu, zxfluxv,
[2469]2634    !   frugs,     agesno,    fsollw,  fsolsw,
2635    !   d_ts,      fevap,     fluxlat, t2m,
2636    !   wfbils,    wfbilo,    fluxt,   fluxu, fluxv,
2637    !
2638    ! Certains ne sont pas utiliser du tout :
2639    !   dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2640    !
[1724]2641
[2469]2642    ! Calcul de l'humidite de saturation au niveau du sol
[1724]2643
2644
[996]2645
[2692]2646    IF (iflag_pbl/=0) THEN
[2240]2647
[2469]2648       !jyg+nrlmd<
[2852]2649!!jyg       IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
2650       IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,10) .ge. 1) THEN
[3179]2651          print *,'debut du splitting de la PBL, wake_s = ', wake_s(:)
2652          print *,'debut du splitting de la PBL, wake_deltat = ', wake_deltat(:,1)
2653          print *,'debut du splitting de la PBL, wake_deltaq = ', wake_deltaq(:,1)
[2469]2654       ENDIF
2655       ! !!
2656       !>jyg+nrlmd
2657       !
2658       !-------gustiness calculation-------!
[3435]2659       !ym : Warning gustiness non inialized for iflag_gusts=2 & iflag_gusts=3
2660       gustiness=0  !ym missing init
2661       
[2469]2662       IF (iflag_gusts==0) THEN
2663          gustiness(1:klon)=0
2664       ELSE IF (iflag_gusts==1) THEN
[3111]2665          gustiness(1:klon)=f_gust_bl*ale_bl(1:klon)+f_gust_wk*ale_wake(1:klon)
2666       ELSE IF (iflag_gusts==2) THEN
2667          gustiness(1:klon)=f_gust_bl*ale_bl_stat(1:klon)+f_gust_wk*ale_wake(1:klon)
[2469]2668          ! ELSE IF (iflag_gusts==2) THEN
2669          !    do i = 1, klon
2670          !       gustiness(i)=f_gust_bl*ale_bl(i)+sigma_wk(i)*f_gust_wk&
2671          !           *ale_wake(i) !! need to make sigma_wk accessible here
2672          !    enddo
2673          ! ELSE IF (iflag_gusts==3) THEN
2674          !    do i = 1, klon
2675          !       gustiness(i)=f_gust_bl*alp_bl(i)+f_gust_wk*alp_wake(i)
2676          !    enddo
2677       ENDIF
[2278]2678
[2469]2679       CALL pbl_surface(  &
[3435]2680            phys_tstep,     date0,     itap,    days_elapsed+1, &
[2469]2681            debut,     lafin, &
2682            longitude_deg, latitude_deg, rugoro,  zrmu0,      &
[4056]2683            sollwdown,    cldt,      &
[3756]2684            rain_fall, snow_fall, solsw,   solswfdiff, sollw,     &
[2469]2685            gustiness,                                &
2686            t_seri,    q_seri,    u_seri,  v_seri,    &
2687                                !nrlmd+jyg<
2688            wake_deltat, wake_deltaq, wake_cstar, wake_s, &
2689                                !>nrlmd+jyg
2690            pplay,     paprs,     pctsrf,             &
2691            ftsol,SFRWL,falb_dir,falb_dif,ustar,u10m,v10m,wstar, &
2692                                !albedo SB <<<
2693            cdragh,    cdragm,  u1,    v1,            &
[3888]2694            beta_aridity, &
[2469]2695                                !albedo SB >>>
2696                                ! albsol1,   albsol2,   sens,    evap,      &
2697            albsol_dir,   albsol_dif,   sens,    evap,   & 
2698                                !albedo SB <<<
2699            albsol3_lic,runoff,   snowhgt,   qsnow, to_ice, sissnow, &
[3817]2700            zxtsol,    zxfluxlat, zt2m,    qsat2m,  zn2mout, &
[2469]2701            d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf, d_t_diss, &
2702                                !nrlmd<
2703                                !jyg<
2704            d_t_vdf_w, d_q_vdf_w, &
2705            d_t_vdf_x, d_q_vdf_x, &
2706            sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
2707                                !>jyg
2708            delta_tsurf,wake_dens, &
2709            cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
2710            kh,kh_x,kh_w, &
2711                                !>nrlmd
2712            coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
2713            slab_wfbils,                 &
2714            qsol,      zq2m,      s_pblh,  s_lcl, &
2715                                !jyg<
2716            s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
2717                                !>jyg
2718            s_capCL,   s_oliqCL,  s_cteiCL,s_pblT, &
2719            s_therm,   s_trmb1,   s_trmb2, s_trmb3, &
2720            zustar, zu10m,     zv10m,   fder, &
[3888]2721            zxqsurf, delta_qsurf,   rh2m,      zxfluxu, zxfluxv, &
[2469]2722            z0m, z0h,     agesno,    fsollw,  fsolsw, &
2723            d_ts,      fevap,     fluxlat, t2m, &
[2670]2724            wfbils, wfbilo, wfevap, wfrain, wfsnow, &
2725            fluxt,   fluxu,  fluxv, &
[2469]2726            dsens,     devap,     zxsnow, &
2727            zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke, &
2728                                !nrlmd+jyg<
[2952]2729            wake_delta_pbl_TKE, &
[2469]2730                                !>nrlmd+jyg
[2952]2731             treedrg )
2732!FC
[2469]2733       !
2734       !  Add turbulent diffusion tendency to the wake difference variables
[2852]2735!!jyg       IF (mod(iflag_pbl_split,2) .NE. 0) THEN
2736       IF (mod(iflag_pbl_split,10) .NE. 0) THEN
[2635]2737!jyg<
2738          d_deltat_vdf(:,:) = d_t_vdf_w(:,:)-d_t_vdf_x(:,:)
2739          d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
2740          CALL add_wake_tend &
[3208]2741             (d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy)
[2635]2742       ELSE
2743          d_deltat_vdf(:,:) = 0.
2744          d_deltaq_vdf(:,:) = 0.
2745!>jyg
[2469]2746       ENDIF
[1624]2747
[2469]2748       !---------------------------------------------------------------------
2749       ! ajout des tendances de la diffusion turbulente
2750       IF (klon_glo==1) THEN
2751          CALL add_pbl_tend &
2752               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
[2799]2753               'vdf',abortphy,flag_inhib_tend,itap)
[2469]2754       ELSE
2755          CALL add_phys_tend &
2756               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
[2812]2757               'vdf',abortphy,flag_inhib_tend,itap,0)
[2469]2758       ENDIF
[3461]2759       CALL prt_enerbil('vdf',itap)
[2469]2760       !--------------------------------------------------------------------
[766]2761
[2692]2762       IF (mydebug) THEN
2763          CALL writefield_phy('u_seri',u_seri,nbp_lev)
2764          CALL writefield_phy('v_seri',v_seri,nbp_lev)
2765          CALL writefield_phy('t_seri',t_seri,nbp_lev)
2766          CALL writefield_phy('q_seri',q_seri,nbp_lev)
2767       ENDIF
[2227]2768
[2469]2769       !albedo SB >>>
2770       albsol1=0.
2771       albsol2=0.
2772       falb1=0.
2773       falb2=0.
[2692]2774       SELECT CASE(nsw)
2775       CASE(2)
[2469]2776          albsol1=albsol_dir(:,1)
2777          albsol2=albsol_dir(:,2)
2778          falb1=falb_dir(:,1,:)
2779          falb2=falb_dir(:,2,:)
[2692]2780       CASE(4)
[2469]2781          albsol1=albsol_dir(:,1)
2782          albsol2=albsol_dir(:,2)*SFRWL(2)+albsol_dir(:,3)*SFRWL(3) &
2783               +albsol_dir(:,4)*SFRWL(4)
2784          albsol2=albsol2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2785          falb1=falb_dir(:,1,:)
2786          falb2=falb_dir(:,2,:)*SFRWL(2)+falb_dir(:,3,:)*SFRWL(3) &
2787               +falb_dir(:,4,:)*SFRWL(4)
2788          falb2=falb2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
[2692]2789       CASE(6)
[2469]2790          albsol1=albsol_dir(:,1)*SFRWL(1)+albsol_dir(:,2)*SFRWL(2) &
2791               +albsol_dir(:,3)*SFRWL(3)
2792          albsol1=albsol1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2793          albsol2=albsol_dir(:,4)*SFRWL(4)+albsol_dir(:,5)*SFRWL(5) &
2794               +albsol_dir(:,6)*SFRWL(6)
2795          albsol2=albsol2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2796          falb1=falb_dir(:,1,:)*SFRWL(1)+falb_dir(:,2,:)*SFRWL(2) &
2797               +falb_dir(:,3,:)*SFRWL(3)
2798          falb1=falb1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2799          falb2=falb_dir(:,4,:)*SFRWL(4)+falb_dir(:,5,:)*SFRWL(5) &
2800               +falb_dir(:,6,:)*SFRWL(6)
2801          falb2=falb2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
[2692]2802       END SELECt
[2469]2803       !albedo SB <<<
[2227]2804
[766]2805
[2469]2806       CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
2807            t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
[1724]2808
[2469]2809    ENDIF
2810    ! =================================================================== c
2811    !   Calcul de Qsat
[881]2812
[2469]2813    DO k = 1, klev
2814       DO i = 1, klon
2815          zx_t = t_seri(i,k)
2816          IF (thermcep) THEN
2817             zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2818             zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2819             zx_qs  = MIN(0.5,zx_qs)
2820             zcor   = 1./(1.-retv*zx_qs)
2821             zx_qs  = zx_qs*zcor
2822          ELSE
2823             !!           IF (zx_t.LT.t_coup) THEN             !jyg
2824             IF (zx_t.LT.rtt) THEN                  !jyg
2825                zx_qs = qsats(zx_t)/pplay(i,k)
2826             ELSE
2827                zx_qs = qsatl(zx_t)/pplay(i,k)
2828             ENDIF
2829          ENDIF
2830          zqsat(i,k)=zx_qs
2831       ENDDO
2832    ENDDO
[959]2833
[2692]2834    IF (prt_level.ge.1) THEN
[2469]2835       write(lunout,*) 'L   qsat (g/kg) avant clouds_gno'
2836       write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
[2692]2837    ENDIF
[2469]2838    !
2839    ! Appeler la convection (au choix)
2840    !
2841    DO k = 1, klev
2842       DO i = 1, klon
2843          conv_q(i,k) = d_q_dyn(i,k)  &
[3435]2844               + d_q_vdf(i,k)/phys_tstep
[2469]2845          conv_t(i,k) = d_t_dyn(i,k)  &
[3435]2846               + d_t_vdf(i,k)/phys_tstep
[2469]2847       ENDDO
2848    ENDDO
2849    IF (check) THEN
2850       za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2851       WRITE(lunout,*) "avantcon=", za
2852    ENDIF
2853    zx_ajustq = .FALSE.
2854    IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2855    IF (zx_ajustq) THEN
2856       DO i = 1, klon
2857          z_avant(i) = 0.0
2858       ENDDO
2859       DO k = 1, klev
2860          DO i = 1, klon
2861             z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
2862                  *(paprs(i,k)-paprs(i,k+1))/RG
2863          ENDDO
2864       ENDDO
2865    ENDIF
[959]2866
[2469]2867    ! Calcule de vitesse verticale a partir de flux de masse verticale
2868    DO k = 1, klev
2869       DO i = 1, klon
2870          omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
[2692]2871       ENDDO
2872    ENDDO
2873
2874    IF (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
[2469]2875         omega(igout, :)
[2707]2876    !
2877    ! Appel de la convection tous les "cvpas"
2878    !
[3150]2879!!jyg    IF (MOD(itapcv,cvpas).EQ.0) THEN
[3151]2880!!    print *,' physiq : itapcv, cvpas, itap-1, cvpas_0 ', &
2881!!                       itapcv, cvpas, itap-1, cvpas_0
[3150]2882    IF (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itap-1,cvpas_0).EQ.0) THEN
[2707]2883
[3134]2884    !
2885    ! Mettre a zero des variables de sortie (pour securite)
2886    !
2887    pmflxr(:,:) = 0.
2888    pmflxs(:,:) = 0.
2889    wdtrainA(:,:) = 0.
[3496]2890    wdtrainS(:,:) = 0.
[3134]2891    wdtrainM(:,:) = 0.
2892    upwd(:,:) = 0.
2893    dnwd(:,:) = 0.
2894    ep(:,:) = 0.
2895    da(:,:)=0.
2896    mp(:,:)=0.
2897    wght_cvfd(:,:)=0.
2898    phi(:,:,:)=0.
2899    phi2(:,:,:)=0.
2900    epmlmMm(:,:,:)=0.
2901    eplaMm(:,:)=0.
2902    d1a(:,:)=0.
2903    dam(:,:)=0.
2904    elij(:,:,:)=0.
2905    ev(:,:)=0.
[3496]2906    qtaa(:,:)=0.
[3134]2907    clw(:,:)=0.
2908    sij(:,:,:)=0.
2909    !
[2469]2910    IF (iflag_con.EQ.1) THEN
2911       abort_message ='reactiver le call conlmd dans physiq.F'
2912       CALL abort_physic (modname,abort_message,1)
[3435]2913       !     CALL conlmd (phys_tstep, paprs, pplay, t_seri, q_seri, conv_q,
[2469]2914       !    .             d_t_con, d_q_con,
2915       !    .             rain_con, snow_con, ibas_con, itop_con)
2916    ELSE IF (iflag_con.EQ.2) THEN
[3435]2917       CALL conflx(phys_tstep, paprs, pplay, t_seri, q_seri, &
[2469]2918            conv_t, conv_q, -evap, omega, &
2919            d_t_con, d_q_con, rain_con, snow_con, &
2920            pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
2921            kcbot, kctop, kdtop, pmflxr, pmflxs)
2922       d_u_con = 0.
2923       d_v_con = 0.
[879]2924
[2469]2925       WHERE (rain_con < 0.) rain_con = 0.
2926       WHERE (snow_con < 0.) snow_con = 0.
2927       DO i = 1, klon
2928          ibas_con(i) = klev+1 - kcbot(i)
2929          itop_con(i) = klev+1 - kctop(i)
2930       ENDDO
2931    ELSE IF (iflag_con.GE.3) THEN
2932       ! nb of tracers for the KE convection:
2933       ! MAF la partie traceurs est faite dans phytrac
2934       ! on met ntra=1 pour limiter les appels mais on peut
2935       ! supprimer les calculs / ftra.
2936       ntra = 1
2937
2938       !=======================================================================
2939       !ajout pour la parametrisation des poches froides: calcul de
[2635]2940       !t_w et t_x: si pas de poches froides, t_w=t_x=t_seri
[2692]2941       IF (iflag_wake>=1) THEN
2942         DO k=1,klev
2943            DO i=1,klon
2944                t_w(i,k) = t_seri(i,k) + (1-wake_s(i))*wake_deltat(i,k)
2945                q_w(i,k) = q_seri(i,k) + (1-wake_s(i))*wake_deltaq(i,k)
2946                t_x(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
2947                q_x(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
2948            ENDDO
2949         ENDDO
2950       ELSE
[4056]2951                t_w(:,:) = t_seri(:,:)
[2635]2952                q_w(:,:) = q_seri(:,:)
2953                t_x(:,:) = t_seri(:,:)
2954                q_x(:,:) = q_seri(:,:)
[2692]2955       ENDIF
[2469]2956       !
2957       !jyg<
2958       ! Perform dry adiabatic adjustment on wake profile
2959       ! The corresponding tendencies are added to the convective tendencies
2960       ! after the call to the convective scheme.
2961       IF (iflag_wake>=1) then
[2882]2962          IF (iflag_adjwk >= 1) THEN
[2469]2963             limbas(:) = 1
[2635]2964             CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
[2309]2965                  d_t_adjwk, d_q_adjwk)
[2638]2966             !
2967             DO k=1,klev
2968                DO i=1,klon
2969                   IF (wake_s(i) .GT. 1.e-3) THEN
2970                      t_w(i,k) = t_w(i,k) + d_t_adjwk(i,k)
2971                      q_w(i,k) = q_w(i,k) + d_q_adjwk(i,k)
2972                      d_deltat_ajs_cv(i,k) = d_t_adjwk(i,k)
2973                      d_deltaq_ajs_cv(i,k) = d_q_adjwk(i,k)
2974                   ELSE
2975                      d_deltat_ajs_cv(i,k) = 0.
2976                      d_deltaq_ajs_cv(i,k) = 0.
2977                   ENDIF
2978                ENDDO
[2469]2979             ENDDO
[2882]2980             IF (iflag_adjwk == 2) THEN
2981               CALL add_wake_tend &
[3208]2982                 (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy)
[2882]2983             ENDIF  ! (iflag_adjwk == 2)
2984          ENDIF  ! (iflag_adjwk >= 1)
[2469]2985       ENDIF ! (iflag_wake>=1)
2986       !>jyg
2987       !
[2638]2988       
2989!!      print *,'physiq. q_w(1,k), q_x(1,k) ', &
2990!!             (k, q_w(1,k), q_x(1,k),k=1,25)
2991
[2513]2992!jyg<
[3435]2993       CALL alpale( debut, itap, phys_tstep, paprs, omega, t_seri,   &
[2513]2994                    alp_offset, it_wape_prescr,  wape_prescr, fip_prescr, &
2995                    ale_bl_prescr, alp_bl_prescr, &
2996                    wake_pe, wake_fip,  &
2997                    Ale_bl, Ale_bl_trig, Alp_bl, &
[2554]2998                    Ale, Alp , Ale_wake, Alp_wake)
[2513]2999!>jyg
3000!
[2469]3001       ! sb, oct02:
3002       ! Schema de convection modularise et vectorise:
3003       ! (driver commun aux versions 3 et 4)
3004       !
3005       IF (ok_cvl) THEN ! new driver for convectL
3006          !
3007          !jyg<
3008          ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3009          ! Calculate the upmost level of deep convection loops: k_upper_cv
3010          !  (near 22 km)
3011          k_upper_cv = klev
[3199]3012          !izero = klon/2+1/klon
3013          !DO k = klev,1,-1
3014          !   IF (pphi(izero,k) > 22.e4) k_upper_cv = k
3015          !ENDDO
3016          ! FH : nouveau calcul base sur un profil global sans quoi
3017          ! le modele etait sensible au decoupage de domaines
[2469]3018          DO k = klev,1,-1
[3199]3019             IF (-7*log(presnivs(k)/presnivs(1)) > 25.) k_upper_cv = k
[2469]3020          ENDDO
3021          IF (prt_level .ge. 5) THEN
3022             Print *, 'upmost level of deep convection loops: k_upper_cv = ', &
3023                  k_upper_cv
3024          ENDIF
3025          !
3026          !>jyg
[4389]3027          IF (type_trac == 'repr') THEN
[2469]3028             nbtr_tmp=ntra
3029          ELSE
3030             nbtr_tmp=nbtr
[2692]3031          ENDIF
[2469]3032          !jyg   iflag_con est dans clesphys
3033          !c          CALL concvl (iflag_con,iflag_clos,
3034          CALL concvl (iflag_clos, &
[3435]3035               phys_tstep, paprs, pplay, k_upper_cv, t_x,q_x, &
[2635]3036               t_w,q_w,wake_s, &
[2469]3037               u_seri,v_seri,tr_seri,nbtr_tmp, &
3038               ALE,ALP, &
3039               sig1,w01, &
3040               d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
3041               rain_con, snow_con, ibas_con, itop_con, sigd, &
[2824]3042               ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
[3496]3043               Ma,mipsh,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
[2469]3044               pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
3045                                ! RomP >>>
3046                                !!     .        pmflxr,pmflxs,da,phi,mp,
3047                                !!     .        ftd,fqd,lalim_conv,wght_th)
[3496]3048               pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,qtaa,clw,elij, &
[2469]3049               ftd,fqd,lalim_conv,wght_th, &
3050               ev, ep,epmlmMm,eplaMm, &
[3496]3051               wdtrainA, wdtrainS, wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
[2481]3052               tau_cld_cv,coefw_cld_cv,epmax_diag)
[2630]3053
[2469]3054          ! RomP <<<
[619]3055
[2469]3056          !IM begin
3057          !       print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
3058          !    .dnwd0(1,1),ftd(1,1),fqd(1,1)
3059          !IM end
3060          !IM cf. FH
3061          clwcon0=qcondc
3062          pmfu(:,:)=upwd(:,:)+dnwd(:,:)
[3150]3063          !
3064          !jyg<
3065          ! If convective tendencies are too large, then call convection
3066          !  every time step
3067          cvpas = cvpas_0
3068          DO k=1,k_upper_cv
3069             DO i=1,klon
[3161]3070               IF (d_t_con(i,k) > 6.721 .AND. d_t_con(i,k) < 6.722 .AND.&
3071                   d_q_con(i,k) > -.0002171 .AND. d_q_con(i,k) < -.0002170) THEN
3072                     dtcon_multistep_max = 3.
3073                     dqcon_multistep_max = 0.02
3074               ENDIF
3075             ENDDO
3076          ENDDO
3077!
3078          DO k=1,k_upper_cv
3079             DO i=1,klon
[3150]3080!!               IF (abs(d_t_con(i,k)) > 0.24 .OR. &
3081!!                   abs(d_q_con(i,k)) > 2.e-2) THEN
3082               IF (abs(d_t_con(i,k)) > dtcon_multistep_max .OR. &
3083                   abs(d_q_con(i,k)) > dqcon_multistep_max) THEN
3084                 cvpas = 1
3085!!                 print *,'physiq1, i,k,d_t_con(i,k),d_q_con(i,k) ', &
3086!!                                   i,k,d_t_con(i,k),d_q_con(i,k)
3087               ENDIF
3088             ENDDO
3089          ENDDO
[3153]3090!!!   Ligne a ne surtout pas remettre sans avoir murement reflechi (jyg)
3091!!!          call bcast(cvpas)
3092!!!   ------------------------------------------------------------
[3150]3093          !>jyg
3094          !
[2692]3095          DO i = 1, klon
[3148]3096             IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+cvpas
[2692]3097          ENDDO
[2469]3098          !
3099          !jyg<
3100          !    Add the tendency due to the dry adjustment of the wake profile
3101          IF (iflag_wake>=1) THEN
[2882]3102            IF (iflag_adjwk == 2) THEN
3103              DO k=1,klev
3104                 DO i=1,klon
[3435]3105                    ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/phys_tstep
3106                    fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/phys_tstep
[2882]3107                    d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
3108                    d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
3109                 ENDDO
3110              ENDDO
3111            ENDIF  ! (iflag_adjwk = 2)
3112          ENDIF   ! (iflag_wake>=1)
[2469]3113          !>jyg
3114          !
3115       ELSE ! ok_cvl
[1412]3116
[2469]3117          ! MAF conema3 ne contient pas les traceurs
[3435]3118          CALL conema3 (phys_tstep, &
[2469]3119               paprs,pplay,t_seri,q_seri, &
3120               u_seri,v_seri,tr_seri,ntra, &
3121               sig1,w01, &
3122               d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
3123               rain_con, snow_con, ibas_con, itop_con, &
3124               upwd,dnwd,dnwd0,bas,top, &
3125               Ma,cape,tvp,rflag, &
3126               pbase &
3127               ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
3128               ,clwcon0)
[524]3129
[2469]3130       ENDIF ! ok_cvl
[524]3131
[2469]3132       !
3133       ! Correction precip
3134       rain_con = rain_con * cvl_corr
3135       snow_con = snow_con * cvl_corr
3136       !
[766]3137
[2469]3138       IF (.NOT. ok_gust) THEN
3139          do i = 1, klon
3140             wd(i)=0.0
3141          enddo
3142       ENDIF
[524]3143
[2469]3144       ! =================================================================== c
3145       ! Calcul des proprietes des nuages convectifs
3146       !
[524]3147
[2469]3148       !   calcul des proprietes des nuages convectifs
3149       clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
3150       IF (iflag_cld_cv == 0) THEN
[2692]3151          CALL clouds_gno &
[2469]3152               (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
3153       ELSE
[2692]3154          CALL clouds_bigauss &
[2469]3155               (klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
3156       ENDIF
[524]3157
[2205]3158
[2469]3159       ! =================================================================== c
[524]3160
[2469]3161       DO i = 1, klon
3162          itop_con(i) = min(max(itop_con(i),1),klev)
3163          ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
3164       ENDDO
[1428]3165
[2469]3166       DO i = 1, klon
[4056]3167          ! C Risi modif: pour éviter pb de dépassement d'indice dans les cas
3168          ! où i n'est pas un point convectif et donc ibas_con(i)=0
3169          ! c'est un pb indépendant des isotopes
3170          if (ibas_con(i) > 0) then
3171             ema_pcb(i)  = paprs(i,ibas_con(i))
3172          else
3173             ema_pcb(i)  = 0.0
3174          endif
[2469]3175       ENDDO
3176       DO i = 1, klon
3177          ! L'idicage de itop_con peut cacher un pb potentiel
3178          ! FH sous la dictee de JYG, CR
3179          ema_pct(i)  = paprs(i,itop_con(i)+1)
[879]3180
[2692]3181          IF (itop_con(i).gt.klev-3) THEN
3182             IF (prt_level >= 9) THEN
[2469]3183                write(lunout,*)'La convection monte trop haut '
3184                write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
[2692]3185             ENDIF
3186          ENDIF
[2469]3187       ENDDO
3188    ELSE IF (iflag_con.eq.0) THEN
3189       write(lunout,*) 'On n appelle pas la convection'
3190       clwcon0=0.
3191       rnebcon0=0.
3192       d_t_con=0.
3193       d_q_con=0.
3194       d_u_con=0.
3195       d_v_con=0.
3196       rain_con=0.
3197       snow_con=0.
3198       bas=1
3199       top=1
3200    ELSE
3201       WRITE(lunout,*) "iflag_con non-prevu", iflag_con
[2692]3202       CALL abort_physic("physiq", "", 1)
[2469]3203    ENDIF
[524]3204
[2469]3205    !     CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
3206    !    .              d_u_con, d_v_con)
[524]3207
[2730]3208!jyg    Reinitialize proba_notrig and itapcv when convection has been called
3209    proba_notrig(:) = 1.
[2707]3210    itapcv = 0
[3150]3211    ENDIF !  (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itapcv,cvpas_0).EQ.0)
[2730]3212!
[2707]3213    itapcv = itapcv+1
[3153]3214    !
3215    ! Compter les steps ou cvpas=1
3216    IF (cvpas == 1) THEN
3217      Ncvpaseq1 = Ncvpaseq1+1
3218    ENDIF
3219    IF (mod(itap,1000) == 0) THEN
3220      print *,' physiq, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
3221    ENDIF
[2707]3222
[2812]3223!!!jyg  Appel diagnostique a add_phys_tend pour tester la conservation de
3224!!!     l'energie dans les courants satures.
3225!!    d_t_con_sat(:,:) = d_t_con(:,:) - ftd(:,:)*dtime
3226!!    d_q_con_sat(:,:) = d_q_con(:,:) - fqd(:,:)*dtime
3227!!    dql_sat(:,:) = (wdtrainA(:,:)+wdtrainM(:,:))*dtime/zmasse(:,:)
3228!!    CALL add_phys_tend(d_u_con, d_v_con, d_t_con_sat, d_q_con_sat, dql_sat,   &
3229!!                     dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,& 
3230!!                     itap, 1)
3231!!    call prt_enerbil('convection_sat',itap)
3232!!
3233!!
[2469]3234    CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
[2812]3235         'convection',abortphy,flag_inhib_tend,itap,0)
[3461]3236    CALL prt_enerbil('convection',itap)
[2235]3237
[2469]3238    !-------------------------------------------------------------------------
[766]3239
[2692]3240    IF (mydebug) THEN
3241       CALL writefield_phy('u_seri',u_seri,nbp_lev)
3242       CALL writefield_phy('v_seri',v_seri,nbp_lev)
3243       CALL writefield_phy('t_seri',t_seri,nbp_lev)
3244       CALL writefield_phy('q_seri',q_seri,nbp_lev)
3245    ENDIF
[766]3246
[2469]3247    IF (check) THEN
3248       za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3249       WRITE(lunout,*)"aprescon=", za
3250       zx_t = 0.0
3251       za = 0.0
3252       DO i = 1, klon
3253          za = za + cell_area(i)/REAL(klon)
3254          zx_t = zx_t + (rain_con(i)+ &
3255               snow_con(i))*cell_area(i)/REAL(klon)
3256       ENDDO
[3435]3257       zx_t = zx_t/za*phys_tstep
[2469]3258       WRITE(lunout,*)"Precip=", zx_t
3259    ENDIF
3260    IF (zx_ajustq) THEN
3261       DO i = 1, klon
3262          z_apres(i) = 0.0
3263       ENDDO
3264       DO k = 1, klev
3265          DO i = 1, klon
3266             z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
3267                  *(paprs(i,k)-paprs(i,k+1))/RG
3268          ENDDO
3269       ENDDO
3270       DO i = 1, klon
[3435]3271          z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*phys_tstep) &
[2469]3272               /z_apres(i)
3273       ENDDO
3274       DO k = 1, klev
3275          DO i = 1, klon
3276             IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
3277                  z_factor(i).LT.(1.0-1.0E-08)) THEN
3278                q_seri(i,k) = q_seri(i,k) * z_factor(i)
3279             ENDIF
3280          ENDDO
3281       ENDDO
3282    ENDIF
3283    zx_ajustq=.FALSE.
[879]3284
[2469]3285    !
3286    !==========================================================================
3287    !RR:Evolution de la poche froide: on ne fait pas de separation wake/env
3288    !pour la couche limite diffuse pour l instant
3289    !
3290    !
3291    ! nrlmd le 22/03/2011---Si on met les poches hors des thermiques
3292    ! il faut rajouter cette tendance calcul\'ee hors des poches
3293    ! froides
3294    !
[2692]3295    IF (iflag_wake>=1) THEN
[2707]3296       !
3297       !
[2730]3298       ! Call wakes every "wkpas" step
3299       !
3300       IF (MOD(itapwk,wkpas).EQ.0) THEN
3301          !
3302          DO k=1,klev
[2469]3303             DO i=1,klon
[2730]3304                dt_dwn(i,k)  = ftd(i,k)
3305                dq_dwn(i,k)  = fqd(i,k)
3306                M_dwn(i,k)   = dnwd0(i,k)
3307                M_up(i,k)    = upwd(i,k)
[3435]3308                dt_a(i,k)    = d_t_con(i,k)/phys_tstep - ftd(i,k)
3309                dq_a(i,k)    = d_q_con(i,k)/phys_tstep - fqd(i,k)
[2469]3310             ENDDO
3311          ENDDO
[2730]3312         
3313          IF (iflag_wake==2) THEN
3314             ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3315             DO k = 1,klev
3316                dt_dwn(:,k)= dt_dwn(:,k)+ &
[3435]3317                     ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/phys_tstep
[2730]3318                dq_dwn(:,k)= dq_dwn(:,k)+ &
[3435]3319                     ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/phys_tstep
[2730]3320             ENDDO
3321          ELSEIF (iflag_wake==3) THEN
3322             ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3323             DO k = 1,klev
3324                DO i=1,klon
3325                   IF (rneb(i,k)==0.) THEN
3326                      ! On ne tient compte des tendances qu'en dehors des
3327                      ! nuages (c'est-\`a-dire a priri dans une region ou
3328                      ! l'eau se reevapore).
3329                      dt_dwn(i,k)= dt_dwn(i,k)+ &
[3435]3330                           ok_wk_lsp(i)*d_t_lsc(i,k)/phys_tstep
[2730]3331                      dq_dwn(i,k)= dq_dwn(i,k)+ &
[3435]3332                           ok_wk_lsp(i)*d_q_lsc(i,k)/phys_tstep
[2730]3333                   ENDIF
3334                ENDDO
3335             ENDDO
3336          ENDIF
3337         
3338          !
3339          !calcul caracteristiques de la poche froide
[3435]3340          CALL calWAKE (iflag_wake_tend, paprs, pplay, phys_tstep, &
[2730]3341               t_seri, q_seri, omega,  &
3342               dt_dwn, dq_dwn, M_dwn, M_up,  &
[3208]3343               dt_a, dq_a, cv_gen,  &
3344               sigd, cin,  &
3345               wake_deltat, wake_deltaq, wake_s, awake_dens, wake_dens,  &
[2730]3346               wake_dth, wake_h,  &
[3000]3347!!               wake_pe, wake_fip, wake_gfl,  &
3348               wake_pe, wake_fip_0, wake_gfl,  &   !! jyg
[2730]3349               d_t_wake, d_q_wake,  &
3350               wake_k, t_x, q_x,  &
3351               wake_omgbdth, wake_dp_omgb,  &
3352               wake_dtKE, wake_dqKE,  &
3353               wake_omg, wake_dp_deltomg,  &
3354               wake_spread, wake_Cstar, d_deltat_wk_gw,  &
[3208]3355               d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk)
[2730]3356          !
3357          !jyg    Reinitialize itapwk when wakes have been called
3358          itapwk = 0
3359       ENDIF !  (MOD(itapwk,wkpas).EQ.0)
[2469]3360       !
[2730]3361       itapwk = itapwk+1
[2469]3362       !
3363       !-----------------------------------------------------------------------
3364       ! ajout des tendances des poches froides
3365       CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
[2812]3366            abortphy,flag_inhib_tend,itap,0)
[3461]3367       CALL prt_enerbil('wake',itap)
[2469]3368       !------------------------------------------------------------------------
[879]3369
[2730]3370       ! Increment Wake state variables
[2635]3371       IF (iflag_wake_tend .GT. 0.) THEN
3372
3373         CALL add_wake_tend &
[3208]3374            (d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk, wake_k, &
[2635]3375             'wake', abortphy)
[3461]3376          CALL prt_enerbil('wake',itap)
[2635]3377       ENDIF   ! (iflag_wake_tend .GT. 0.)
[3179]3378       !
3379       IF (prt_level .GE. 10) THEN
3380         print *,' physiq, after calwake, wake_s: ',wake_s(:)
3381         print *,' physiq, after calwake, wake_deltat: ',wake_deltat(:,1)
3382         print *,' physiq, after calwake, wake_deltaq: ',wake_deltaq(:,1)
3383       ENDIF
[2635]3384
[3000]3385       IF (iflag_alp_wk_cond .GT. 0.) THEN
3386
[3435]3387         CALL alpale_wk(phys_tstep, cell_area, wake_k, wake_s, wake_dens, wake_fip_0, &
[3000]3388                        wake_fip)
3389       ELSE
3390         wake_fip(:) = wake_fip_0(:)
3391       ENDIF   ! (iflag_alp_wk_cond .GT. 0.)
3392
[2692]3393    ENDIF  ! (iflag_wake>=1)
[2469]3394    !
3395    !===================================================================
3396    ! Convection seche (thermiques ou ajustement)
3397    !===================================================================
3398    !
[2692]3399    CALL stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
[2469]3400         ,seuil_inversion,weak_inversion,dthmin)
[878]3401
3402
3403
[2469]3404    d_t_ajsb(:,:)=0.
3405    d_q_ajsb(:,:)=0.
3406    d_t_ajs(:,:)=0.
3407    d_u_ajs(:,:)=0.
3408    d_v_ajs(:,:)=0.
3409    d_q_ajs(:,:)=0.
3410    clwcon0th(:,:)=0.
3411    !
3412    !      fm_therm(:,:)=0.
3413    !      entr_therm(:,:)=0.
3414    !      detr_therm(:,:)=0.
3415    !
[2692]3416    IF (prt_level>9) WRITE(lunout,*) &
[2469]3417         'AVANT LA CONVECTION SECHE , iflag_thermals=' &
3418         ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
[2692]3419    IF (iflag_thermals<0) THEN
[2469]3420       !  Rien
3421       !  ====
[2692]3422       IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
[541]3423
[878]3424
[2692]3425    ELSE
[878]3426
[2469]3427       !  Thermiques
3428       !  ==========
[2692]3429       IF (prt_level>9) WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
[2469]3430            ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
[878]3431
3432
[2469]3433       !cc nrlmd le 10/04/2012
3434       DO k=1,klev+1
3435          DO i=1,klon
3436             pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
3437             pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
3438             pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
3439             pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
[2159]3440          ENDDO
[2469]3441       ENDDO
3442       !cc fin nrlmd le 10/04/2012
[1403]3443
[2692]3444       IF (iflag_thermals>=1) THEN
[2469]3445          !jyg<
[2852]3446!!       IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3447       IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
[2469]3448             !  Appel des thermiques avec les profils exterieurs aux poches
3449             DO k=1,klev
3450                DO i=1,klon
3451                   t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
3452                   q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
[2606]3453                   u_therm(i,k) = u_seri(i,k)
3454                   v_therm(i,k) = v_seri(i,k)
[2469]3455                ENDDO
3456             ENDDO
3457          ELSE
3458             !  Appel des thermiques avec les profils moyens
3459             DO k=1,klev
3460                DO i=1,klon
3461                   t_therm(i,k) = t_seri(i,k)
3462                   q_therm(i,k) = q_seri(i,k)
[2606]3463                   u_therm(i,k) = u_seri(i,k)
3464                   v_therm(i,k) = v_seri(i,k)
[2469]3465                ENDDO
3466             ENDDO
3467          ENDIF
3468          !>jyg
[2692]3469          CALL calltherm(pdtphys &
[2469]3470               ,pplay,paprs,pphi,weak_inversion &
[2606]3471                        ! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
3472               ,u_therm,v_therm,t_therm,q_therm,zqsat,debut &  !jyg
[2469]3473               ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
3474               ,fm_therm,entr_therm,detr_therm &
3475               ,zqasc,clwcon0th,lmax_th,ratqscth &
3476               ,ratqsdiff,zqsatth &
3477                                !on rajoute ale et alp, et les
3478                                !caracteristiques de la couche alim
3479               ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
3480               ,ztv,zpspsk,ztla,zthl &
3481                                !cc nrlmd le 10/04/2012
3482               ,pbl_tke_input,pctsrf,omega,cell_area &
3483               ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
3484               ,n2,s2,ale_bl_stat &
3485               ,therm_tke_max,env_tke_max &
3486               ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
3487               ,alp_bl_conv,alp_bl_stat &
3488                                !cc fin nrlmd le 10/04/2012
3489               ,zqla,ztva )
3490          !
3491          !jyg<
[2852]3492!!jyg          IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3493          IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
[2469]3494             !  Si les thermiques ne sont presents que hors des
3495             !  poches, la tendance moyenne associ\'ee doit etre
3496             !  multipliee par la fraction surfacique qu'ils couvrent.
3497             DO k=1,klev
3498                DO i=1,klon
3499                   !
[2635]3500                   d_deltat_the(i,k) = - d_t_ajs(i,k)
3501                   d_deltaq_the(i,k) = - d_q_ajs(i,k)
[2469]3502                   !
3503                   d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
3504                   d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
3505                   d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
3506                   d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
3507                   !
3508                ENDDO
3509             ENDDO
[2606]3510          !
[3180]3511             IF (ok_bug_split_th) THEN
3512               CALL add_wake_tend &
[3208]3513                   (d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy)
[3180]3514             ELSE
3515               CALL add_wake_tend &
[3208]3516                   (d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wake_k, 'the', abortphy)
[3180]3517             ENDIF
[3461]3518             CALL prt_enerbil('the',itap)
[2638]3519          !
[2852]3520          ENDIF  ! (mod(iflag_pbl_split/10,10) .GE. 1)
[2638]3521          !
[2606]3522          CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,  &
[2812]3523                             dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
[3461]3524          CALL prt_enerbil('thermals',itap)
[2606]3525          !
[2513]3526!
[3435]3527          CALL alpale_th( phys_tstep, lmax_th, t_seri, cell_area,  &
[2513]3528                          cin, s2, n2,  &
3529                          ale_bl_trig, ale_bl_stat, ale_bl,  &
[2556]3530                          alp_bl, alp_bl_stat, &
[3208]3531                          proba_notrig, random_notrig, cv_gen)
[2635]3532          !>jyg
[1638]3533
[2554]3534          ! ------------------------------------------------------------------
3535          ! Transport de la TKE par les panaches thermiques.
3536          ! FH : 2010/02/01
3537          !     if (iflag_pbl.eq.10) then
3538          !     call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
3539          !    s           rg,paprs,pbl_tke)
3540          !     endif
3541          ! -------------------------------------------------------------------
3542
[2692]3543          DO i=1,klon
[2469]3544             !           zmax_th(i)=pphi(i,lmax_th(i))/rg
3545             !CR:04/05/12:correction calcul zmax
3546             zmax_th(i)=zmax0(i)
[2692]3547          ENDDO
[1507]3548
[2692]3549       ENDIF
[878]3550
[2469]3551       !  Ajustement sec
3552       !  ==============
[878]3553
[2469]3554       ! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
3555       ! a partir du sommet des thermiques.
3556       ! Dans le cas contraire, on demarre au niveau 1.
[878]3557
[2692]3558       IF (iflag_thermals>=13.or.iflag_thermals<=0) THEN
[878]3559
[2692]3560          IF (iflag_thermals.eq.0) THEN
3561             IF (prt_level>9) WRITE(lunout,*)'ajsec'
[2469]3562             limbas(:)=1
[2692]3563          ELSE
[2469]3564             limbas(:)=lmax_th(:)
[2692]3565          ENDIF
[878]3566
[2469]3567          ! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
3568          ! pour des test de convergence numerique.
3569          ! Le nouveau ajsec est a priori mieux, meme pour le cas
3570          ! iflag_thermals = 0 (l'ancienne version peut faire des tendances
3571          ! non nulles numeriquement pour des mailles non concernees.
[878]3572
[2692]3573          IF (iflag_thermals==0) THEN
[2469]3574             ! Calling adjustment alone (but not the thermal plume model)
3575             CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
3576                  , d_t_ajsb, d_q_ajsb)
[2692]3577          ELSE IF (iflag_thermals>0) THEN
[2469]3578             ! Calling adjustment above the top of thermal plumes
3579             CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
3580                  , d_t_ajsb, d_q_ajsb)
[2692]3581          ENDIF
[878]3582
[2469]3583          !--------------------------------------------------------------------
3584          ! ajout des tendances de l'ajustement sec ou des thermiques
3585          CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
[2812]3586               'ajsb',abortphy,flag_inhib_tend,itap,0)
[3461]3587          CALL prt_enerbil('ajsb',itap)
[2469]3588          d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
3589          d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
[904]3590
[2469]3591          !---------------------------------------------------------------------
[878]3592
[2692]3593       ENDIF
[524]3594
[2692]3595    ENDIF
[2469]3596    !
3597    !===================================================================
3598    ! Computation of ratqs, the width (normalized) of the subrid scale
3599    ! water distribution
[4009]3600
3601    tke_dissip_ave(:,:)=0.
3602    l_mix_ave(:,:)=0.
3603    wprime_ave(:,:)=0.
3604
3605    DO nsrf = 1, nbsrf
3606       DO i = 1, klon
3607          tke_dissip_ave(i,:) = tke_dissip_ave(i,:) + tke_dissip(i,:,nsrf)*pctsrf(i,nsrf)
3608          l_mix_ave(i,:) = l_mix_ave(i,:) + l_mix(i,:,nsrf)*pctsrf(i,nsrf)
3609          wprime_ave(i,:) = wprime_ave(i,:) + wprime(i,:,nsrf)*pctsrf(i,nsrf)
3610       ENDDO
3611    ENDDO
3612
[2469]3613    CALL  calcratqs(klon,klev,prt_level,lunout,        &
3614         iflag_ratqs,iflag_con,iflag_cld_th,pdtphys,  &
[2534]3615         ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
[3856]3616         tau_ratqs,fact_cldcon,wake_s, wake_deltaq,   &
[2469]3617         ptconv,ptconvth,clwcon0th, rnebcon0th,     &
[4009]3618         paprs,pplay,t_seri,q_seri, qtc_cv, sigt_cv, zqsat, &
3619         pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave,t2m,q2m,fm_therm, &
[3856]3620         ratqs,ratqsc,ratqs_inter)
[1032]3621
[2469]3622    !
3623    ! Appeler le processus de condensation a grande echelle
3624    ! et le processus de precipitation
3625    !-------------------------------------------------------------------------
3626    IF (prt_level .GE.10) THEN
3627       print *,'itap, ->fisrtilp ',itap
3628    ENDIF
[4056]3629    !
[3999]3630
3631    picefra(:,:)=0.
3632
3633    IF (ok_new_lscp) THEN
3634
[4062]3635    !--mise à jour de flight_m et flight_h2o dans leur module
3636    IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
3637      CALL airplane(debut,pphis,pplay,paprs,t_seri)
3638    ENDIF
[4059]3639
[4380]3640    CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay, &
[3999]3641         t_seri, q_seri,ptconv,ratqs, &
[4380]3642         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, rneb_seri, &
[4412]3643         radocond, picefra, rain_lsc, snow_lsc, &
[3999]3644         frac_impa, frac_nucl, beta_prec_fisrt, &
3645         prfl, psfl, rhcl,  &
3646         zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
[4380]3647         iflag_ice_thermo, ok_ice_sursat, zqsatl, zqsats, &
3648         qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
3649         Tcontr, qcontr, qcontr2, fcontrN, fcontrP )
[3999]3650
3651    ELSE
[4059]3652
[3435]3653    CALL fisrtilp(phys_tstep,paprs,pplay, &
[2469]3654         t_seri, q_seri,ptconv,ratqs, &
[4412]3655         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, radocond, &
[2469]3656         rain_lsc, snow_lsc, &
3657         pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3658         frac_impa, frac_nucl, beta_prec_fisrt, &
3659         prfl, psfl, rhcl,  &
3660         zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3661         iflag_ice_thermo)
[4059]3662
[3999]3663    ENDIF
[4056]3664    !
[2469]3665    WHERE (rain_lsc < 0) rain_lsc = 0.
3666    WHERE (snow_lsc < 0) snow_lsc = 0.
[766]3667
[2799]3668!+JLD
3669!    write(*,9000) 'phys lsc',"enerbil: bil_q, bil_e,",rain_lsc+snow_lsc &
3670!        & ,((rcw-rcpd)*rain_lsc + (rcs-rcpd)*snow_lsc)*t_seri(1,1)-rlvtt*rain_lsc+rlstt*snow_lsc &
3671!        & ,rain_lsc,snow_lsc
3672!    write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
3673!-JLD
[2469]3674    CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
[2812]3675         'lsc',abortphy,flag_inhib_tend,itap,0)
[3461]3676    CALL prt_enerbil('lsc',itap)
[2613]3677    rain_num(:)=0.
[2657]3678    DO k = 1, klev
[2613]3679       DO i = 1, klon
3680          IF (ql_seri(i,k)>oliqmax) THEN
3681             rain_num(i)=rain_num(i)+(ql_seri(i,k)-oliqmax)*zmasse(i,k)/pdtphys
3682             ql_seri(i,k)=oliqmax
3683          ENDIF
3684       ENDDO
3685    ENDDO
[4098]3686    IF (nqo >= 3) THEN
[2657]3687    DO k = 1, klev
3688       DO i = 1, klon
3689          IF (qs_seri(i,k)>oicemax) THEN
3690             rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
3691             qs_seri(i,k)=oicemax
3692          ENDIF
3693       ENDDO
3694    ENDDO
3695    ENDIF
[2613]3696
[2524]3697    !---------------------------------------------------------------------------
[2469]3698    DO k = 1, klev
3699       DO i = 1, klon
3700          cldfra(i,k) = rneb(i,k)
3701          !CR: a quoi ca sert? Faut-il ajouter qs_seri?
[4412]3702          IF (.NOT.new_oliq) radocond(i,k) = ql_seri(i,k)
[2469]3703       ENDDO
3704    ENDDO
3705    IF (check) THEN
3706       za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3707       WRITE(lunout,*)"apresilp=", za
3708       zx_t = 0.0
3709       za = 0.0
3710       DO i = 1, klon
3711          za = za + cell_area(i)/REAL(klon)
3712          zx_t = zx_t + (rain_lsc(i) &
3713               + snow_lsc(i))*cell_area(i)/REAL(klon)
3714       ENDDO
[3435]3715       zx_t = zx_t/za*phys_tstep
[2469]3716       WRITE(lunout,*)"Precip=", zx_t
3717    ENDIF
[766]3718
[2692]3719    IF (mydebug) THEN
3720       CALL writefield_phy('u_seri',u_seri,nbp_lev)
3721       CALL writefield_phy('v_seri',v_seri,nbp_lev)
3722       CALL writefield_phy('t_seri',t_seri,nbp_lev)
3723       CALL writefield_phy('q_seri',q_seri,nbp_lev)
3724    ENDIF
[524]3725
[2469]3726    !
3727    !-------------------------------------------------------------------
3728    !  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3729    !-------------------------------------------------------------------
[524]3730
[2469]3731    ! 1. NUAGES CONVECTIFS
3732    !
3733    !IM cf FH
3734    !     IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
3735    IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
3736       snow_tiedtke=0.
3737       !     print*,'avant calcul de la pseudo precip '
3738       !     print*,'iflag_cld_th',iflag_cld_th
[2692]3739       IF (iflag_cld_th.eq.-1) THEN
[2469]3740          rain_tiedtke=rain_con
[2692]3741       ELSE
[2469]3742          !       print*,'calcul de la pseudo precip '
3743          rain_tiedtke=0.
3744          !         print*,'calcul de la pseudo precip 0'
[2692]3745          DO k=1,klev
3746             DO i=1,klon
3747                IF (d_q_con(i,k).lt.0.) THEN
[2469]3748                   rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
3749                        *(paprs(i,k)-paprs(i,k+1))/rg
[2692]3750                ENDIF
3751             ENDDO
3752          ENDDO
3753       ENDIF
[2469]3754       !
3755       !     call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3756       !
[524]3757
[2469]3758       ! Nuages diagnostiques pour Tiedtke
3759       CALL diagcld1(paprs,pplay, &
3760                                !IM cf FH. rain_con,snow_con,ibas_con,itop_con,
3761            rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
3762            diafra,dialiq)
3763       DO k = 1, klev
3764          DO i = 1, klon
3765             IF (diafra(i,k).GT.cldfra(i,k)) THEN
[4412]3766                radocond(i,k) = dialiq(i,k)
[2469]3767                cldfra(i,k) = diafra(i,k)
3768             ENDIF
3769          ENDDO
3770       ENDDO
[524]3771
[2469]3772    ELSE IF (iflag_cld_th.ge.3) THEN
3773       !  On prend pour les nuages convectifs le max du calcul de la
3774       !  convection et du calcul du pas de temps precedent diminue d'un facteur
3775       !  facttemps
3776       facteur = pdtphys *facttemps
[2692]3777       DO k=1,klev
3778          DO i=1,klon
[2469]3779             rnebcon(i,k)=rnebcon(i,k)*facteur
[2692]3780             IF (rnebcon0(i,k)*clwcon0(i,k).GT.rnebcon(i,k)*clwcon(i,k)) THEN
[2469]3781                rnebcon(i,k)=rnebcon0(i,k)
3782                clwcon(i,k)=clwcon0(i,k)
[2692]3783             ENDIF
3784          ENDDO
3785       ENDDO
[2469]3786
3787       !   On prend la somme des fractions nuageuses et des contenus en eau
[524]3788
[2692]3789       IF (iflag_cld_th>=5) THEN
[1411]3790
[2692]3791          DO k=1,klev
[2469]3792             ptconvth(:,k)=fm_therm(:,k+1)>0.
[2692]3793          ENDDO
[1496]3794
[2692]3795          IF (iflag_coupl==4) THEN
[1496]3796
[2469]3797             ! Dans le cas iflag_coupl==4, on prend la somme des convertures
3798             ! convectives et lsc dans la partie des thermiques
3799             ! Le controle par iflag_coupl est peut etre provisoire.
[2692]3800             DO k=1,klev
3801                DO i=1,klon
3802                   IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
[4412]3803                      radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
[2469]3804                      cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
[2692]3805                   ELSE IF (ptconv(i,k)) THEN
[2469]3806                      cldfra(i,k)=rnebcon(i,k)
[4412]3807                      radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
[2692]3808                   ENDIF
3809                ENDDO
3810             ENDDO
[1496]3811
[2692]3812          ELSE IF (iflag_coupl==5) THEN
3813             DO k=1,klev
3814                DO i=1,klon
[2469]3815                   cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
[4412]3816                   radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
[2692]3817                ENDDO
3818             ENDDO
[1525]3819
[2692]3820          ELSE
[1525]3821
[2469]3822             ! Si on est sur un point touche par la convection
3823             ! profonde et pas par les thermiques, on prend la
3824             ! couverture nuageuse et l'eau nuageuse de la convection
3825             ! profonde.
[1411]3826
[2469]3827             !IM/FH: 2011/02/23
3828             ! definition des points sur lesquels ls thermiques sont actifs
[1496]3829
[2692]3830             DO k=1,klev
3831                DO i=1,klon
3832                   IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
[2469]3833                      cldfra(i,k)=rnebcon(i,k)
[4412]3834                      radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
[2692]3835                   ENDIF
3836                ENDDO
3837             ENDDO
[1496]3838
[2692]3839          ENDIF
[1496]3840
[2692]3841       ELSE
[1496]3842
[2469]3843          ! Ancienne version
3844          cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
[4412]3845          radocond(:,:)=radocond(:,:)+rnebcon(:,:)*clwcon(:,:)
[2692]3846       ENDIF
[1411]3847
[2469]3848    ENDIF
[1507]3849
[2469]3850    !     plulsc(:)=0.
3851    !     do k=1,klev,-1
3852    !        do i=1,klon
3853    !              zzz=prfl(:,k)+psfl(:,k)
3854    !           if (.not.ptconvth.zzz.gt.0.)
3855    !        enddo prfl, psfl,
3856    !     enddo
3857    !
3858    ! 2. NUAGES STARTIFORMES
3859    !
3860    IF (ok_stratus) THEN
3861       CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3862       DO k = 1, klev
3863          DO i = 1, klon
3864             IF (diafra(i,k).GT.cldfra(i,k)) THEN
[4412]3865                radocond(i,k) = dialiq(i,k)
[2469]3866                cldfra(i,k) = diafra(i,k)
3867             ENDIF
3868          ENDDO
3869       ENDDO
3870    ENDIF
3871    !
3872    ! Precipitation totale
3873    !
3874    DO i = 1, klon
3875       rain_fall(i) = rain_con(i) + rain_lsc(i)
3876       snow_fall(i) = snow_con(i) + snow_lsc(i)
3877    ENDDO
3878    !
3879    ! Calculer l'humidite relative pour diagnostique
3880    !
3881    DO k = 1, klev
3882       DO i = 1, klon
3883          zx_t = t_seri(i,k)
3884          IF (thermcep) THEN
3885             !!           if (iflag_ice_thermo.eq.0) then                 !jyg
3886             zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3887             !!           else                                            !jyg
3888             !!           zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t))      !jyg
3889             !!           endif                                           !jyg
3890             zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3891             zx_qs  = MIN(0.5,zx_qs)
3892             zcor   = 1./(1.-retv*zx_qs)
3893             zx_qs  = zx_qs*zcor
3894          ELSE
3895             !!           IF (zx_t.LT.t_coup) THEN             !jyg
3896             IF (zx_t.LT.rtt) THEN                  !jyg
3897                zx_qs = qsats(zx_t)/pplay(i,k)
3898             ELSE
3899                zx_qs = qsatl(zx_t)/pplay(i,k)
3900             ENDIF
3901          ENDIF
3902          zx_rh(i,k) = q_seri(i,k)/zx_qs
[3784]3903            IF (iflag_ice_thermo .GT. 0) THEN
[3780]3904          zx_rhl(i,k) = q_seri(i,k)/(qsatl(zx_t)/pplay(i,k))
3905          zx_rhi(i,k) = q_seri(i,k)/(qsats(zx_t)/pplay(i,k))
[3784]3906            ENDIF
[2469]3907          zqsat(i,k)=zx_qs
3908       ENDDO
3909    ENDDO
[782]3910
[2469]3911    !IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3912    !   equivalente a 2m (tpote) pour diagnostique
3913    !
3914    DO i = 1, klon
3915       tpot(i)=zt2m(i)*(100000./paprs(i,1))**RKAPPA
3916       IF (thermcep) THEN
3917          IF(zt2m(i).LT.RTT) then
3918             Lheat=RLSTT
3919          ELSE
3920             Lheat=RLVTT
3921          ENDIF
3922       ELSE
3923          IF (zt2m(i).LT.RTT) THEN
3924             Lheat=RLSTT
3925          ELSE
3926             Lheat=RLVTT
3927          ENDIF
3928       ENDIF
3929       tpote(i) = tpot(i)*      &
3930            EXP((Lheat *qsat2m(i))/(RCPD*zt2m(i)))
3931    ENDDO
[524]3932
[4389]3933    IF (ANY(type_trac == ['inca','inco'])) THEN ! ModThL
[524]3934#ifdef INCA
[2469]3935       CALL VTe(VTphysiq)
3936       CALL VTb(VTinca)
3937       calday = REAL(days_elapsed + 1) + jH_cur
[524]3938
[3435]3939       CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap)
[3563]3940       CALL AEROSOL_METEO_CALC( &
3941            calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3942            prfl,psfl,pctsrf,cell_area, &
3943            latitude_deg,longitude_deg,u10m,v10m)
[524]3944
[2469]3945       zxsnow_dummy(:) = 0.0
[625]3946
[2469]3947       CALL chemhook_begin (calday, &
3948            days_elapsed+1, &
3949            jH_cur, &
3950            pctsrf(1,1), &
3951            latitude_deg, &
3952            longitude_deg, &
3953            cell_area, &
3954            paprs, &
3955            pplay, &
3956            coefh(1:klon,1:klev,is_ave), &
3957            pphi, &
3958            t_seri, &
3959            u, &
3960            v, &
[3613]3961            rot, &
[2469]3962            wo(:, :, 1), &
3963            q_seri, &
3964            zxtsol, &
[3613]3965            zt2m, &
[2469]3966            zxsnow_dummy, &
3967            solsw, &
3968            albsol1, &
3969            rain_fall, &
3970            snow_fall, &
3971            itop_con, &
3972            ibas_con, &
3973            cldfra, &
3974            nbp_lon, &
3975            nbp_lat-1, &
[3872]3976            tr_seri(:,:,1+nqCO2:nbtr), &
[2469]3977            ftsol, &
3978            paprs, &
3979            cdragh, &
3980            cdragm, &
3981            pctsrf, &
3982            pdtphys, &
3983            itap)
[616]3984
[2469]3985       CALL VTe(VTinca)
3986       CALL VTb(VTphysiq)
[3865]3987#endif
3988    ENDIF !type_trac = inca or inco
[4389]3989    IF (type_trac == 'repr') THEN
[3666]3990#ifdef REPROBUS
3991    !CALL chemtime_rep(itap+itau_phy-1, date0, dtime, itap)
3992    CALL chemtime_rep(itap+itau_phy-1, date0, phys_tstep, itap)
3993#endif
3994    ENDIF
[2618]3995
[2469]3996    !
[2618]3997    ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3998    !
3999    IF (MOD(itaprad,radpas).EQ.0) THEN
[959]4000
[2618]4001       !
4002       !jq - introduce the aerosol direct and first indirect radiative forcings
4003       !jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
[2738]4004       IF (flag_aerosol .GT. 0) THEN
[2618]4005          IF (iflag_rrtm .EQ. 0) THEN !--old radiation
4006             IF (.NOT. aerosol_couple) THEN
4007                !
4008                CALL readaerosol_optic( &
[3630]4009                     debut, flag_aerosol, itap, jD_cur-jD_ref, &
[2618]4010                     pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
4011                     mass_solu_aero, mass_solu_aero_pi,  &
4012                     tau_aero, piz_aero, cg_aero,  &
4013                     tausum_aero, tau3d_aero)
4014             ENDIF
4015          ELSE                       ! RRTM radiation
4016             IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
4017                abort_message='config_inca=aero et rrtm=1 impossible'
[2692]4018                CALL abort_physic(modname,abort_message,1)
[2618]4019             ELSE
4020                !
4021#ifdef CPP_RRTM
4022                IF (NSW.EQ.6) THEN
[2738]4023                   !--new aerosol properties SW and LW
[2618]4024                   !
[2753]4025#ifdef CPP_Dust
4026                   !--SPL aerosol model
4027                   CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, &
4028                        tr_seri, mass_solu_aero, mass_solu_aero_pi,  &
4029                        tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,  &
4030                        tausum_aero, tau3d_aero)
4031#else
4032                   !--climatologies or INCA aerosols
[3480]4033                   CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, &
[3630]4034                        flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
[2618]4035                        pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
4036                        tr_seri, mass_solu_aero, mass_solu_aero_pi,  &
4037                        tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,  &
[2854]4038                        tausum_aero, drytausum_aero, tau3d_aero)
[2753]4039#endif
[3274]4040
4041                   IF (flag_aerosol .EQ. 7) THEN
4042                      CALL MACv2SP(pphis,pplay,paprs,longitude_deg,latitude_deg,  &
4043                                   tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
4044                   ENDIF
4045
[2738]4046                   !
[2618]4047                ELSE IF (NSW.EQ.2) THEN
4048                   !--for now we use the old aerosol properties
4049                   !
4050                   CALL readaerosol_optic( &
[3630]4051                        debut, flag_aerosol, itap, jD_cur-jD_ref, &
[2618]4052                        pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
4053                        mass_solu_aero, mass_solu_aero_pi,  &
4054                        tau_aero, piz_aero, cg_aero,  &
4055                        tausum_aero, tau3d_aero)
4056                   !
4057                   !--natural aerosols
4058                   tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
4059                   piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
4060                   cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
4061                   !--all aerosols
4062                   tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
4063                   piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
4064                   cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
[2738]4065                   !
4066                   !--no LW optics
4067                   tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
4068                   !
[2618]4069                ELSE
4070                   abort_message='Only NSW=2 or 6 are possible with ' &
4071                        // 'aerosols and iflag_rrtm=1'
[2692]4072                   CALL abort_physic(modname,abort_message,1)
[2618]4073                ENDIF
4074#else
4075                abort_message='You should compile with -rrtm if running ' &
4076                     // 'with iflag_rrtm=1'
[2692]4077                CALL abort_physic(modname,abort_message,1)
[2618]4078#endif
4079                !
4080             ENDIF
4081          ENDIF
[2738]4082       ELSE   !--flag_aerosol = 0
[2618]4083          tausum_aero(:,:,:) = 0.
[2854]4084          drytausum_aero(:,:) = 0.
[2640]4085          mass_solu_aero(:,:) = 0.
4086          mass_solu_aero_pi(:,:) = 0.
[2618]4087          IF (iflag_rrtm .EQ. 0) THEN !--old radiation
4088             tau_aero(:,:,:,:) = 1.e-15
4089             piz_aero(:,:,:,:) = 1.
4090             cg_aero(:,:,:,:)  = 0.
4091          ELSE
4092             tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
4093             tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
4094             piz_aero_sw_rrtm(:,:,:,:) = 1.0
4095             cg_aero_sw_rrtm(:,:,:,:)  = 0.0
4096          ENDIF
4097       ENDIF
4098       !
[2994]4099       !--WMO criterion to determine tropopause
[3123]4100       CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
[2994]4101       !
[2618]4102       !--STRAT AEROSOL
4103       !--updates tausum_aero,tau_aero,piz_aero,cg_aero
4104       IF (flag_aerosol_strat.GT.0) THEN
4105          IF (prt_level .GE.10) THEN
4106             PRINT *,'appel a readaerosolstrat', mth_cur
4107          ENDIF
4108          IF (iflag_rrtm.EQ.0) THEN
4109           IF (flag_aerosol_strat.EQ.1) THEN
4110             CALL readaerosolstrato(debut)
4111           ELSE
4112             abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
4113             CALL abort_physic(modname,abort_message,1)
4114           ENDIF
4115          ELSE
[2009]4116#ifdef CPP_RRTM
[2690]4117#ifndef CPP_StratAer
4118          !--prescribed strat aerosols
4119          !--only in the case of non-interactive strat aerosols
[2618]4120            IF (flag_aerosol_strat.EQ.1) THEN
4121             CALL readaerosolstrato1_rrtm(debut)
4122            ELSEIF (flag_aerosol_strat.EQ.2) THEN
[3480]4123             CALL readaerosolstrato2_rrtm(debut, ok_volcan)
[2618]4124            ELSE
4125             abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
4126             CALL abort_physic(modname,abort_message,1)
4127            ENDIF
[2690]4128#endif
[2618]4129#else
4130             abort_message='You should compile with -rrtm if running ' &
4131                  // 'with iflag_rrtm=1'
4132             CALL abort_physic(modname,abort_message,1)
4133#endif
4134          ENDIF
[3567]4135       ELSE
4136          tausum_aero(:,:,id_STRAT_phy) = 0.
[2618]4137       ENDIF
[2690]4138!
4139#ifdef CPP_RRTM
4140#ifdef CPP_StratAer
[2692]4141       !--compute stratospheric mask
[3123]4142       CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
[2690]4143       !--interactive strat aerosols
4144       CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
4145#endif
4146#endif
[2618]4147       !--fin STRAT AEROSOL
4148       !     
4149
4150       ! Calculer les parametres optiques des nuages et quelques
4151       ! parametres pour diagnostiques:
4152       !
4153       IF (aerosol_couple.AND.config_inca=='aero') THEN
4154          mass_solu_aero(:,:)    = ccm(:,:,1)
4155          mass_solu_aero_pi(:,:) = ccm(:,:,2)
[2692]4156       ENDIF
[2618]4157
4158       IF (ok_newmicro) then
[3908]4159! AI          IF (iflag_rrtm.NE.0) THEN
4160          IF (iflag_rrtm.EQ.1) THEN
[2618]4161#ifdef CPP_RRTM
4162             IF (ok_cdnc.AND.NRADLP.NE.3) THEN
[2469]4163             abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
4164                  // 'pour ok_cdnc'
[2618]4165             CALL abort_physic(modname,abort_message,1)
4166             ENDIF
[2009]4167#else
4168
[2618]4169             abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
4170             CALL abort_physic(modname,abort_message,1)
[2009]4171#endif
[2618]4172          ENDIF
[3274]4173          CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
[4412]4174               paprs, pplay, t_seri, radocond, picefra, cldfra, &
[2618]4175               cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
4176               flwp, fiwp, flwc, fiwc, &
4177               mass_solu_aero, mass_solu_aero_pi, &
[4114]4178               cldtaupi, latitude_deg, re, fl, ref_liq, ref_ice, &
[2618]4179               ref_liq_pi, ref_ice_pi)
4180       ELSE
4181          CALL nuage (paprs, pplay, &
[4412]4182               t_seri, radocond, picefra, cldfra, cldtau, cldemi, &
[2618]4183               cldh, cldl, cldm, cldt, cldq, &
4184               ok_aie, &
4185               mass_solu_aero, mass_solu_aero_pi, &
4186               bl95_b0, bl95_b1, &
4187               cldtaupi, re, fl)
[2469]4188       ENDIF
4189       !
[2618]4190       !IM betaCRF
[2469]4191       !
[2618]4192       cldtaurad   = cldtau
4193       cldtaupirad = cldtaupi
4194       cldemirad   = cldemi
4195       cldfrarad   = cldfra
4196
[2469]4197       !
[2618]4198       IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
4199           lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
4200          !
4201          ! global
4202          !
[3048]4203!IM 251017 begin
[3317]4204!               print*,'physiq betaCRF global zdtime=',zdtime
[3048]4205!IM 251017 end
[2618]4206          DO k=1, klev
4207             DO i=1, klon
4208                IF (pplay(i,k).GE.pfree) THEN
[2469]4209                   beta(i,k) = beta_pbl
[2618]4210                ELSE
[2469]4211                   beta(i,k) = beta_free
[2618]4212                ENDIF
4213                IF (mskocean_beta) THEN
[2469]4214                   beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
[2618]4215                ENDIF
[2469]4216                cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
4217                cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4218                cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
4219                cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
[2618]4220             ENDDO
4221          ENDDO
4222          !
4223       ELSE
4224          !
4225          ! regional
4226          !
4227          DO k=1, klev
4228             DO i=1,klon
4229                !
4230                IF (longitude_deg(i).ge.lon1_beta.AND. &
4231                    longitude_deg(i).le.lon2_beta.AND. &
4232                    latitude_deg(i).le.lat1_beta.AND.  &
4233                    latitude_deg(i).ge.lat2_beta) THEN
4234                   IF (pplay(i,k).GE.pfree) THEN
4235                      beta(i,k) = beta_pbl
4236                   ELSE
4237                      beta(i,k) = beta_free
4238                   ENDIF
4239                   IF (mskocean_beta) THEN
4240                      beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
4241                   ENDIF
4242                   cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
4243                   cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4244                   cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
4245                   cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
4246                ENDIF
[2469]4247             !
[2618]4248             ENDDO
[2469]4249          ENDDO
4250       !
[2618]4251       ENDIF
[766]4252
[2618]4253       !lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
4254       IF (ok_chlorophyll) THEN
[2469]4255          print*,"-- reading chlorophyll"
[2618]4256          CALL readchlorophyll(debut)
4257       ENDIF
[1863]4258
[2524]4259!--if ok_suntime_rrtm we use ancillay data for RSUN
4260!--previous values are therefore overwritten
4261!--this is needed for CMIP6 runs
4262!--and only possible for new radiation scheme
4263       IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
[2525]4264#ifdef CPP_RRTM
[2524]4265         CALL read_rsun_rrtm(debut)
[2525]4266#endif
[2524]4267       ENDIF
4268
[2692]4269       IF (mydebug) THEN
4270          CALL writefield_phy('u_seri',u_seri,nbp_lev)
4271          CALL writefield_phy('v_seri',v_seri,nbp_lev)
4272          CALL writefield_phy('t_seri',t_seri,nbp_lev)
4273          CALL writefield_phy('q_seri',q_seri,nbp_lev)
4274       ENDIF
[2524]4275
[2469]4276       !
4277       !sonia : If Iflag_radia >=2, pertubation of some variables
4278       !input to radiation (DICE)
4279       !
4280       IF (iflag_radia .ge. 2) THEN
4281          zsav_tsol (:) = zxtsol(:)
[2692]4282          CALL perturb_radlwsw(zxtsol,iflag_radia)
[2469]4283       ENDIF
[2328]4284
[2469]4285       IF (aerosol_couple.AND.config_inca=='aero') THEN
[959]4286#ifdef INCA
[2469]4287          CALL radlwsw_inca  &
[3338]4288               (chemistry_couple, kdlon,kflev,dist, rmu0, fract, solaire, &
[2469]4289               paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
[2684]4290               size(wo,3), wo, &
[2469]4291               cldfrarad, cldemirad, cldtaurad, &
4292               heat,heat0,cool,cool0,albpla, &
4293               topsw,toplw,solsw,sollw, &
4294               sollwdown, &
4295               topsw0,toplw0,solsw0,sollw0, &
4296               lwdn0, lwdn, lwup0, lwup,  &
4297               swdn0, swdn, swup0, swup, &
4298               ok_ade, ok_aie, &
4299               tau_aero, piz_aero, cg_aero, &
4300               topswad_aero, solswad_aero, &
4301               topswad0_aero, solswad0_aero, &
4302               topsw_aero, topsw0_aero, &
4303               solsw_aero, solsw0_aero, &
4304               cldtaupirad, &
4305               topswai_aero, solswai_aero)
[955]4306#endif
[2469]4307       ELSE
4308          !
4309          !IM calcul radiatif pour le cas actuel
4310          !
4311          RCO2 = RCO2_act
4312          RCH4 = RCH4_act
4313          RN2O = RN2O_act
4314          RCFC11 = RCFC11_act
4315          RCFC12 = RCFC12_act
[3450]4316          !
4317          !--interactive CO2 in ppm from carbon cycle
[4298]4318          IF (carbon_cycle_rad) RCO2=RCO2_glo
[2469]4319          !
4320          IF (prt_level .GE.10) THEN
4321             print *,' ->radlwsw, number 1 '
4322          ENDIF
4323          !
4324          CALL radlwsw &
4325               (dist, rmu0, fract,  &
4326                                !albedo SB >>>
4327                                !      paprs, pplay,zxtsol,albsol1, albsol2,  &
4328               paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif,  &
4329                                !albedo SB <<<
4330               t_seri,q_seri,wo, &
4331               cldfrarad, cldemirad, cldtaurad, &
[3989]4332               ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie,  ok_volcan, flag_volc_surfstrat, &
4333               flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
[2469]4334               tau_aero, piz_aero, cg_aero, &
4335               tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4336               ! Rajoute par OB pour RRTM
4337               tau_aero_lw_rrtm, &
[3630]4338               cldtaupirad, &
[3048]4339!              zqsat, flwcrad, fiwcrad, &
[2469]4340               zqsat, flwc, fiwc, &
4341               ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4342               heat,heat0,cool,cool0,albpla, &
[3479]4343               heat_volc,cool_volc, &
[3756]4344               topsw,toplw,solsw,solswfdiff,sollw, &
[2469]4345               sollwdown, &
4346               topsw0,toplw0,solsw0,sollw0, &
[3106]4347               lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup,  &
[3082]4348               swdnc0, swdn0, swdn, swupc0, swup0, swup, &
[2469]4349               topswad_aero, solswad_aero, &
4350               topswai_aero, solswai_aero, &
4351               topswad0_aero, solswad0_aero, &
4352               topsw_aero, topsw0_aero, &
4353               solsw_aero, solsw0_aero, &
4354               topswcf_aero, solswcf_aero, &
4355                                !-C. Kleinschmitt for LW diagnostics
4356               toplwad_aero, sollwad_aero,&
4357               toplwai_aero, sollwai_aero, &
4358               toplwad0_aero, sollwad0_aero,&
4359                                !-end
4360               ZLWFT0_i, ZFLDN0, ZFLUP0, &
[3117]4361               ZSWFT0_i, ZFSDN0, ZFSUP0)
[879]4362
[3048]4363          !lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
4364          !schemes
4365          toplw = toplw + betalwoff * (toplw0 - toplw)
4366          sollw = sollw + betalwoff * (sollw0 - sollw)
4367          lwdn = lwdn + betalwoff * (lwdn0 - lwdn)
4368          lwup = lwup + betalwoff * (lwup0 - lwup)
4369          sollwdown(:)= sollwdown(:) + betalwoff *(-1.*ZFLDN0(:,1) - &
4370                        sollwdown(:))
4371          cool = cool + betalwoff * (cool0 - cool)
4372 
[2679]4373#ifndef CPP_XIOS
[4056]4374          !
[2469]4375          !IM 2eme calcul radiatif pour le cas perturbe ou au moins un
4376          !IM des taux doit etre different du taux actuel
4377          !IM Par defaut on a les taux perturbes egaux aux taux actuels
4378          !
[2989]4379          IF (RCO2_per.NE.RCO2_act.OR. &
4380              RCH4_per.NE.RCH4_act.OR. &
4381              RN2O_per.NE.RN2O_act.OR. &
4382              RCFC11_per.NE.RCFC11_act.OR. &
4383              RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
[3149]4384#endif
[2989]4385   !
[2692]4386          IF (ok_4xCO2atm) THEN
[2469]4387                !
4388                RCO2 = RCO2_per
4389                RCH4 = RCH4_per
4390                RN2O = RN2O_per
4391                RCFC11 = RCFC11_per
4392                RCFC12 = RCFC12_per
4393                !
4394                IF (prt_level .GE.10) THEN
4395                   print *,' ->radlwsw, number 2 '
4396                ENDIF
4397                !
4398                CALL radlwsw &
4399                     (dist, rmu0, fract,  &
4400                                !albedo SB >>>
4401                                !      paprs, pplay,zxtsol,albsol1, albsol2,  &
4402                     paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4403                                !albedo SB <<<
4404                     t_seri,q_seri,wo, &
[2640]4405                     cldfrarad, cldemirad, cldtaurad, &
[3989]4406                     ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie,  ok_volcan, flag_volc_surfstrat, &
4407                     flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
[2469]4408                     tau_aero, piz_aero, cg_aero, &
4409                     tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4410                                ! Rajoute par OB pour RRTM
4411                     tau_aero_lw_rrtm, &
[3630]4412                     cldtaupi, &
[3048]4413!                    zqsat, flwcrad, fiwcrad, &
[2469]4414                     zqsat, flwc, fiwc, &
4415                     ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4416                     heatp,heat0p,coolp,cool0p,albplap, &
[3479]4417                     heat_volc,cool_volc, &
[3756]4418                     topswp,toplwp,solswp,solswfdiffp,sollwp, &
[2469]4419                     sollwdownp, &
4420                     topsw0p,toplw0p,solsw0p,sollw0p, &
[3106]4421                     lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp,  &
[3082]4422                     swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp, &
[2469]4423                     topswad_aerop, solswad_aerop, &
4424                     topswai_aerop, solswai_aerop, &
4425                     topswad0_aerop, solswad0_aerop, &
4426                     topsw_aerop, topsw0_aerop, &
4427                     solsw_aerop, solsw0_aerop, &
4428                     topswcf_aerop, solswcf_aerop, &
4429                                !-C. Kleinschmitt for LW diagnostics
4430                     toplwad_aerop, sollwad_aerop,&
4431                     toplwai_aerop, sollwai_aerop, &
4432                     toplwad0_aerop, sollwad0_aerop,&
4433                                !-end
4434                     ZLWFT0_i, ZFLDN0, ZFLUP0, &
[3117]4435                     ZSWFT0_i, ZFSDN0, ZFSUP0)
[3987]4436          ENDIF !ok_4xCO2atm
[2469]4437       ENDIF ! aerosol_couple
4438       itaprad = 0
4439       !
4440       !  If Iflag_radia >=2, reset pertubed variables
4441       !
4442       IF (iflag_radia .ge. 2) THEN
4443          zxtsol(:) = zsav_tsol (:)
4444       ENDIF
4445    ENDIF ! MOD(itaprad,radpas)
4446    itaprad = itaprad + 1
[879]4447
[2469]4448    IF (iflag_radia.eq.0) THEN
4449       IF (prt_level.ge.9) THEN
4450          PRINT *,'--------------------------------------------------'
4451          PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
4452          PRINT *,'>>>>           heat et cool mis a zero '
4453          PRINT *,'--------------------------------------------------'
[2692]4454       ENDIF
[2469]4455       heat=0.
4456       cool=0.
4457       sollw=0.   ! MPL 01032011
4458       solsw=0.
4459       radsol=0.
4460       swup=0.    ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
4461       swup0=0.
4462       lwup=0.
4463       lwup0=0.
4464       lwdn=0.
4465       lwdn0=0.
[2692]4466    ENDIF
[782]4467
[2469]4468    !
4469    ! Calculer radsol a l'exterieur de radlwsw
4470    ! pour prendre en compte le cycle diurne
4471    ! recode par Olivier Boucher en sept 2015
4472    !
4473    radsol=solsw*swradcorr+sollw
[2618]4474
[2692]4475    IF (ok_4xCO2atm) THEN
[2469]4476       radsolp=solswp*swradcorr+sollwp
[2692]4477    ENDIF
[2359]4478
[2469]4479    !
4480    ! Ajouter la tendance des rayonnements (tous les pas)
4481    ! avec une correction pour le cycle diurne dans le SW
4482    !
[2359]4483
[2469]4484    DO k=1, klev
[3435]4485       d_t_swr(:,k)=swradcorr(:)*heat(:,k)*phys_tstep/RDAY
4486       d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)*phys_tstep/RDAY
4487       d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY
4488       d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY
[2469]4489    ENDDO
[2194]4490
[2812]4491    CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
[3461]4492    CALL prt_enerbil('SW',itap)
[2812]4493    CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
[3461]4494    CALL prt_enerbil('LW',itap)
[1863]4495
[2469]4496    !
[2692]4497    IF (mydebug) THEN
4498       CALL writefield_phy('u_seri',u_seri,nbp_lev)
4499       CALL writefield_phy('v_seri',v_seri,nbp_lev)
4500       CALL writefield_phy('t_seri',t_seri,nbp_lev)
4501       CALL writefield_phy('q_seri',q_seri,nbp_lev)
4502    ENDIF
[1863]4503
[2469]4504    ! Calculer l'hydrologie de la surface
4505    !
4506    !      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
4507    !     .            agesno, ftsol,fqsurf,fsnow, ruis)
4508    !
[1001]4509
[2469]4510    !
4511    ! Calculer le bilan du sol et la derive de temperature (couplage)
4512    !
4513    DO i = 1, klon
4514       !         bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
4515       ! a la demande de JLD
4516       bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
4517    ENDDO
4518    !
4519    !moddeblott(jan95)
4520    ! Appeler le programme de parametrisation de l'orographie
4521    ! a l'echelle sous-maille:
4522    !
4523    IF (prt_level .GE.10) THEN
4524       print *,' call orography ? ', ok_orodr
4525    ENDIF
4526    !
4527    IF (ok_orodr) THEN
4528       !
4529       !  selection des points pour lesquels le shema est actif:
4530       igwd=0
4531       DO i=1,klon
4532          itest(i)=0
[4458]4533          zrel_oro(i)=zstd(i)/(max(zsig(i),1.E-8)*sqrt(cell_area(i)))
4534          !zrel_oro: relative mountain height wrt relief explained by mean slope
4535          ! -> condition on zrel_oro can deactivate the drag on tilted planar terrains
[4352]4536          !    such as ice sheets (work by V. Wiener)
[4458]4537          ! zpmm_orodr_t and zstd_orodr_t are activation thresholds set by F. Lott to
4538          ! earn computation time but they are not physical.
4539          IF (((zpic(i)-zmea(i)).GT.zpmm_orodr_t).AND.(zstd(i).GT.zstd_orodr_t).AND.(zrel_oro(i).GE.zrel_oro_t)) THEN
[2469]4540             itest(i)=1
4541             igwd=igwd+1
4542             idx(igwd)=i
4543          ENDIF
4544       ENDDO
4545       !        igwdim=MAX(1,igwd)
4546       !
4547       IF (ok_strato) THEN
[1863]4548
[3435]4549          CALL drag_noro_strato(0,klon,klev,phys_tstep,paprs,pplay, &
[2469]4550               zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4551               igwd,idx,itest, &
4552               t_seri, u_seri, v_seri, &
4553               zulow, zvlow, zustrdr, zvstrdr, &
4554               d_t_oro, d_u_oro, d_v_oro)
[1863]4555
[2469]4556       ELSE
[3435]4557          CALL drag_noro(klon,klev,phys_tstep,paprs,pplay, &
[2469]4558               zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4559               igwd,idx,itest, &
4560               t_seri, u_seri, v_seri, &
4561               zulow, zvlow, zustrdr, zvstrdr, &
4562               d_t_oro, d_u_oro, d_v_oro)
4563       ENDIF
4564       !
4565       !  ajout des tendances
4566       !-----------------------------------------------------------------------
4567       ! ajout des tendances de la trainee de l'orographie
4568       CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
[2812]4569            abortphy,flag_inhib_tend,itap,0)
[3461]4570       CALL prt_enerbil('oro',itap)
[2469]4571       !----------------------------------------------------------------------
4572       !
4573    ENDIF ! fin de test sur ok_orodr
4574    !
[2692]4575    IF (mydebug) THEN
4576       CALL writefield_phy('u_seri',u_seri,nbp_lev)
4577       CALL writefield_phy('v_seri',v_seri,nbp_lev)
4578       CALL writefield_phy('t_seri',t_seri,nbp_lev)
4579       CALL writefield_phy('q_seri',q_seri,nbp_lev)
4580    ENDIF
[1001]4581
[2469]4582    IF (ok_orolf) THEN
4583       !
4584       !  selection des points pour lesquels le shema est actif:
4585       igwd=0
4586       DO i=1,klon
4587          itest(i)=0
[4458]4588          !zrel_oro: relative mountain height wrt relief explained by mean slope
4589          ! -> condition on zrel_oro can deactivate the lifting on tilted planar terrains
[4352]4590          !    such as ice sheets (work by V. Wiener)
[4458]4591          zrel_oro(i)=zstd(i)/(max(zsig(i),1.E-8)*sqrt(cell_area(i)))
4592          IF (((zpic(i)-zmea(i)).GT.zpmm_orolf_t).AND.(zrel_oro(i).GE.zrel_oro_t)) THEN
[2469]4593             itest(i)=1
4594             igwd=igwd+1
4595             idx(igwd)=i
4596          ENDIF
4597       ENDDO
4598       !        igwdim=MAX(1,igwd)
4599       !
4600       IF (ok_strato) THEN
[1001]4601
[3435]4602          CALL lift_noro_strato(klon,klev,phys_tstep,paprs,pplay, &
[2469]4603               latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
4604               igwd,idx,itest, &
4605               t_seri, u_seri, v_seri, &
4606               zulow, zvlow, zustrli, zvstrli, &
4607               d_t_lif, d_u_lif, d_v_lif               )
[2333]4608
[2469]4609       ELSE
[3435]4610          CALL lift_noro(klon,klev,phys_tstep,paprs,pplay, &
[2469]4611               latitude_deg,zmea,zstd,zpic, &
4612               itest, &
4613               t_seri, u_seri, v_seri, &
4614               zulow, zvlow, zustrli, zvstrli, &
4615               d_t_lif, d_u_lif, d_v_lif)
4616       ENDIF
[1638]4617
[2469]4618       ! ajout des tendances de la portance de l'orographie
4619       CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
[2812]4620            'lif', abortphy,flag_inhib_tend,itap,0)
[3461]4621       CALL prt_enerbil('lif',itap)
[2469]4622    ENDIF ! fin de test sur ok_orolf
[1638]4623
[2469]4624    IF (ok_hines) then
4625       !  HINES GWD PARAMETRIZATION
4626       east_gwstress=0.
4627       west_gwstress=0.
4628       du_gwd_hines=0.
4629       dv_gwd_hines=0.
[3435]4630       CALL hines_gwd(klon, klev, phys_tstep, paprs, pplay, latitude_deg, t_seri, &
[2469]4631            u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
4632            du_gwd_hines, dv_gwd_hines)
4633       zustr_gwd_hines=0.
4634       zvstr_gwd_hines=0.
4635       DO k = 1, klev
[3435]4636          zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/phys_tstep &
[2469]4637               * (paprs(:, k)-paprs(:, k+1))/rg
[3435]4638          zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/phys_tstep &
[2469]4639               * (paprs(:, k)-paprs(:, k+1))/rg
4640       ENDDO
[1001]4641
[2469]4642       d_t_hin(:, :)=0.
4643       CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
[2812]4644            dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
[3461]4645       CALL prt_enerbil('hin',itap)
[2469]4646    ENDIF
[2333]4647
[2469]4648    IF (.not. ok_hines .and. ok_gwd_rando) then
[3435]4649       ! ym missing init for east_gwstress & west_gwstress -> added in phys_local_var_mod
4650       CALL acama_GWD_rando(PHYS_TSTEP, pplay, latitude_deg, t_seri, u_seri, &
[2469]4651            v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
4652            dv_gwd_front, east_gwstress, west_gwstress)
4653       zustr_gwd_front=0.
4654       zvstr_gwd_front=0.
4655       DO k = 1, klev
[3435]4656          zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/phys_tstep &
[2469]4657               * (paprs(:, k)-paprs(:, k+1))/rg
[3435]4658          zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/phys_tstep &
[2469]4659               * (paprs(:, k)-paprs(:, k+1))/rg
4660       ENDDO
[644]4661
[2469]4662       CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
[2812]4663            paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
[3461]4664       CALL prt_enerbil('front_gwd_rando',itap)
[2469]4665    ENDIF
[1938]4666
[2692]4667    IF (ok_gwd_rando) THEN
[3435]4668       CALL FLOTT_GWD_rando(PHYS_TSTEP, pplay, t_seri, u_seri, v_seri, &
[2469]4669            rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
4670            du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
4671       CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
[2812]4672            paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
[3461]4673       CALL prt_enerbil('flott_gwd_rando',itap)
[2469]4674       zustr_gwd_rando=0.
4675       zvstr_gwd_rando=0.
4676       DO k = 1, klev
[3435]4677          zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/phys_tstep &
[2469]4678               * (paprs(:, k)-paprs(:, k+1))/rg
[3435]4679          zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/phys_tstep &
[2469]4680               * (paprs(:, k)-paprs(:, k+1))/rg
4681       ENDDO
[2692]4682    ENDIF
[766]4683
[2469]4684    ! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
[1279]4685
[2692]4686    IF (mydebug) THEN
4687       CALL writefield_phy('u_seri',u_seri,nbp_lev)
4688       CALL writefield_phy('v_seri',v_seri,nbp_lev)
4689       CALL writefield_phy('t_seri',t_seri,nbp_lev)
4690       CALL writefield_phy('q_seri',q_seri,nbp_lev)
4691    ENDIF
[2136]4692
[2469]4693    DO i = 1, klon
4694       zustrph(i)=0.
4695       zvstrph(i)=0.
4696    ENDDO
4697    DO k = 1, klev
4698       DO i = 1, klon
[3435]4699          zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/phys_tstep* &
[2469]4700               (paprs(i,k)-paprs(i,k+1))/rg
[3435]4701          zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/phys_tstep* &
[2469]4702               (paprs(i,k)-paprs(i,k+1))/rg
4703       ENDDO
4704    ENDDO
4705    !
4706    !IM calcul composantes axiales du moment angulaire et couple des montagnes
4707    !
4708    IF (is_sequential .and. ok_orodr) THEN
4709       CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
4710            ra,rg,romega, &
4711            latitude_deg,longitude_deg,pphis, &
4712            zustrdr,zustrli,zustrph, &
4713            zvstrdr,zvstrli,zvstrph, &
4714            paprs,u,v, &
4715            aam, torsfc)
4716    ENDIF
4717    !IM cf. FLott END
4718    !DC Calcul de la tendance due au methane
[3461]4719    IF (ok_qch4) THEN
[2469]4720       CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
4721       ! ajout de la tendance d'humidite due au methane
[3435]4722       d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*phys_tstep
[2801]4723       CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
[2812]4724            'q_ch4', abortphy,flag_inhib_tend,itap,0)
[3435]4725       d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/phys_tstep
[2692]4726    ENDIF
[2469]4727    !
4728    !
[2897]4729
4730!===============================================================
4731!            Additional tendency of TKE due to orography
4732!===============================================================
4733!
4734! Inititialization
4735!------------------
4736
4737       addtkeoro=0   
4738       CALL getin_p('addtkeoro',addtkeoro)
4739     
4740       IF (prt_level.ge.5) &
4741            print*,'addtkeoro', addtkeoro
4742           
4743       alphatkeoro=1.   
4744       CALL getin_p('alphatkeoro',alphatkeoro)
4745       alphatkeoro=min(max(0.,alphatkeoro),1.)
4746
[3461]4747       smallscales_tkeoro=.FALSE.   
[2897]4748       CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
4749
4750
[3461]4751       dtadd(:,:)=0.
4752       duadd(:,:)=0.
4753       dvadd(:,:)=0.
[2897]4754
4755! Choices for addtkeoro:
4756!      ** 0 no TKE tendency from orography   
4757!      ** 1 we include a fraction alphatkeoro of the whole tendency duoro
4758!      ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
4759!
4760
4761       IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
4762!      -------------------------------------------
4763
4764
4765       !  selection des points pour lesquels le schema est actif:
4766
4767
4768  IF (addtkeoro .EQ. 1 ) THEN
4769
4770            duadd(:,:)=alphatkeoro*d_u_oro(:,:)
4771            dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
4772
4773  ELSE IF (addtkeoro .EQ. 2) THEN
4774
[3461]4775     IF (smallscales_tkeoro) THEN
[2897]4776       igwd=0
4777       DO i=1,klon
4778          itest(i)=0
4779! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
4780! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
4781! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
[4458]4782          IF ((zstd(i).GT.1.0) .AND.(zrel_oro(i).GE.zrel_oro_t)) THEN
[2897]4783             itest(i)=1
4784             igwd=igwd+1
4785             idx(igwd)=i
4786          ENDIF
4787       ENDDO
4788
4789     ELSE
4790
4791       igwd=0
4792       DO i=1,klon
4793          itest(i)=0
[4458]4794        IF (((zpic(i)-zmea(i)).GT.zpmm_orodr_t).AND.(zstd(i).GT.zstd_orodr_t).AND.(zrel_oro(i).GE.zrel_oro_t)) THEN
[2897]4795             itest(i)=1
4796             igwd=igwd+1
4797             idx(igwd)=i
[3461]4798        ENDIF
[2897]4799       ENDDO
4800
[3461]4801     ENDIF
[2897]4802
[3461]4803     CALL drag_noro_strato(addtkeoro,klon,klev,phys_tstep,paprs,pplay, &
[2897]4804               zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4805               igwd,idx,itest, &
4806               t_seri, u_seri, v_seri, &
4807               zulow, zvlow, zustrdr, zvstrdr, &
4808               d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
4809
[3461]4810     zustrdr(:)=0.
4811     zvstrdr(:)=0.
4812     zulow(:)=0.
4813     zvlow(:)=0.
[2897]4814
[3461]4815     duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
4816     dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
4817  ENDIF
[2897]4818
4819
4820   ! TKE update from subgrid temperature and wind tendencies
4821   !----------------------------------------------------------
4822    forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4823
4824
[3198]4825    CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
[3888]4826   !
4827   ! Prevent pbl_tke_w from becoming negative
4828    wake_delta_pbl_tke(:,:,:) = max(wake_delta_pbl_tke(:,:,:), -pbl_tke(:,:,:))
4829   !
[2897]4830
4831       ENDIF
4832!      -----
4833!===============================================================
4834
4835
[2469]4836    !====================================================================
4837    ! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4838    !====================================================================
4839    ! Abderrahmane 24.08.09
4840
4841    IF (ok_cosp) THEN
4842       ! adeclarer
[1279]4843#ifdef CPP_COSP
[3435]4844       IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
[1279]4845
[2469]4846          IF (prt_level .GE.10) THEN
4847             print*,'freq_cosp',freq_cosp
4848          ENDIF
4849          mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4850          !       print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4851          !     s        ref_liq,ref_ice
[3435]4852          CALL phys_cosp(itap,phys_tstep,freq_cosp, &
[2469]4853               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
[2794]4854               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
[2469]4855               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4856               JrNt,ref_liq,ref_ice, &
4857               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4858               zu10m,zv10m,pphis, &
4859               zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4860               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4861               prfl(:,1:klev),psfl(:,1:klev), &
4862               pmflxr(:,1:klev),pmflxs(:,1:klev), &
4863               mr_ozone,cldtau, cldemi)
[1412]4864
[2469]4865          !     L         calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4866          !     L          cfaddbze,clcalipso2,dbze,cltlidarradar,
4867          !     M          clMISR,
4868          !     R          clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4869          !     I          tauisccp,albisccp,meantbisccp,meantbclrisccp)
[1279]4870
[2469]4871       ENDIF
[3370]4872#endif
[1279]4873
[3370]4874#ifdef CPP_COSP2
[3435]4875       IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
[3370]4876
4877          IF (prt_level .GE.10) THEN
4878             print*,'freq_cosp',freq_cosp
4879          ENDIF
4880          mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4881                 print*,'Dans physiq.F avant appel '
4882          !     s        ref_liq,ref_ice
[3435]4883          CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
[3370]4884               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4885               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4886               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4887               JrNt,ref_liq,ref_ice, &
4888               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4889               zu10m,zv10m,pphis, &
4890               zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4891               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4892               prfl(:,1:klev),psfl(:,1:klev), &
4893               pmflxr(:,1:klev),pmflxs(:,1:klev), &
4894               mr_ozone,cldtau, cldemi)
4895       ENDIF
[1279]4896#endif
[3370]4897
[3491]4898#ifdef CPP_COSPV2
4899       IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
[3726]4900!        IF (MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
[3491]4901
4902          IF (prt_level .GE.10) THEN
4903             print*,'freq_cosp',freq_cosp
4904          ENDIF
[3726]4905           DO k = 1, klev
4906             DO i = 1, klon
4907               phicosp(i,k) = pphi(i,k) + pphis(i)
4908             ENDDO
4909           ENDDO
[3491]4910          mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4911                 print*,'Dans physiq.F avant appel '
4912          !     s        ref_liq,ref_ice
4913          CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
4914               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4915               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4916               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4917               JrNt,ref_liq,ref_ice, &
4918               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4919               zu10m,zv10m,pphis, &
4920               zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4921               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4922               prfl(:,1:klev),psfl(:,1:klev), &
4923               pmflxr(:,1:klev),pmflxs(:,1:klev), &
4924               mr_ozone,cldtau, cldemi)
4925       ENDIF
4926#endif
4927
[2469]4928    ENDIF  !ok_cosp
[2580]4929
4930
4931! Marine
4932
4933  IF (ok_airs) then
4934
[3435]4935  IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/phys_tstep)).EQ.0) THEN
[2692]4936     write(*,*) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
4937     CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
4938        & map_prop_hc,map_prop_hist,&
4939        & map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
4940        & map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
4941        & map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
4942        & map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
4943        & map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
4944        & map_ntot,map_hc,map_hist,&
4945        & map_Cb,map_ThCi,map_Anv,&
4946        & alt_tropo )
[2580]4947  ENDIF
4948
4949  ENDIF  ! ok_airs
4950
4951
[2469]4952    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4953    !AA
4954    !AA Installation de l'interface online-offline pour traceurs
4955    !AA
4956    !====================================================================
4957    !   Calcul  des tendances traceurs
4958    !====================================================================
4959    !
[959]4960
[4389]4961    IF (type_trac == 'repr') THEN
[3666]4962!MM pas d'impact, car on recupere q_seri,tr_seri,t_seri via phys_local_var_mod
4963!MM                               dans Reprobus
[2469]4964       sh_in(:,:) = q_seri(:,:)
[3666]4965#ifdef REPROBUS
4966       d_q_rep(:,:) = 0.
4967       d_ql_rep(:,:) = 0.
4968       d_qi_rep(:,:) = 0.
4969#endif
[2469]4970    ELSE
4971       sh_in(:,:) = qx(:,:,ivap)
[4098]4972       IF (nqo >= 3) THEN
[3861]4973          ch_in(:,:) = qx(:,:,iliq) + qx(:,:,isol)
4974       ELSE
4975          ch_in(:,:) = qx(:,:,iliq)
4976       ENDIF
[2692]4977    ENDIF
[1565]4978
[2630]4979#ifdef CPP_Dust
[3776]4980    !  Avec SPLA, iflag_phytrac est forcé =1
4981    CALL       phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con,       &  ! I
[2630]4982                      pdtphys,ftsol,                                   &  ! I
4983                      t,q_seri,paprs,pplay,RHcl,                  &  ! I
4984                      pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,          &  ! I
4985                      coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1,                 &  ! I
4986                      u_seri, v_seri, latitude_deg, longitude_deg,  &
4987                      pphis,pctsrf,pmflxr,pmflxs,prfl,psfl,            &  ! I
4988                      da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij,     &  ! I
4989                      epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con,      &  ! I
4990                      ev,wdtrainA,  wdtrainM,wght_cvfd,              &  ! I
4991                      fm_therm, entr_therm, rneb,                      &  ! I
4992                      beta_prec_fisrt,beta_prec, & !I
4993                      zu10m,zv10m,wstar,ale_bl,ale_wake,               &  ! I
4994                      d_tr_dyn,tr_seri)
4995
4996#else
[3776]4997    IF (iflag_phytrac == 1 ) THEN
4998      CALL phytrac ( &
[2469]4999         itap,     days_elapsed+1,    jH_cur,   debut, &
[3435]5000         lafin,    phys_tstep,     u, v,     t, &
[2469]5001         paprs,    pplay,     pmfu,     pmfd, &
5002         pen_u,    pde_u,     pen_d,    pde_d, &
5003         cdragh,   coefh(1:klon,1:klev,is_ave),   fm_therm, entr_therm, &
5004         u1,       v1,        ftsol,    pctsrf, &
5005         zustar,   zu10m,     zv10m, &
5006         wstar(:,is_ave),    ale_bl,         ale_wake, &
5007         latitude_deg, longitude_deg, &
5008         frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
5009         presnivs, pphis,     pphi,     albsol1, &
[2784]5010         sh_in,   ch_in,    rhcl,      cldfra,   rneb, &
[4412]5011         diafra,   radocond,    itop_con, ibas_con, &
[2469]5012         pmflxr,   pmflxs,    prfl,     psfl, &
5013         da,       phi,       mp,       upwd, &
5014         phi2,     d1a,       dam,      sij, wght_cvfd, &        !<<RomP+RL
5015         wdtrainA, wdtrainM,  sigd,     clw,elij, &   !<<RomP
5016         ev,       ep,        epmlmMm,  eplaMm, &     !<<RomP
5017         dnwd,     aerosol_couple,      flxmass_w, &
5018         tau_aero, piz_aero,  cg_aero,  ccm, &
5019         rfname, &
5020         d_tr_dyn, &                                 !<<RomP
[3418]5021         tr_seri, init_source)
[3666]5022#ifdef REPROBUS
5023
5024
5025          print*,'avt add phys rep',abortphy
5026
5027     CALL add_phys_tend &
5028            (du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,paprs,&
5029             'rep',abortphy,flag_inhib_tend,itap,0)
5030        IF (abortphy==1) Print*,'ERROR ABORT REP'
5031
5032          print*,'apr add phys rep',abortphy
5033
[2630]5034#endif
[3776]5035    ENDIF    ! (iflag_phytrac=1)
[3666]5036
5037#endif
[3776]5038    !ENDIF    ! (iflag_phytrac=1)
[524]5039
[2469]5040    IF (offline) THEN
[524]5041
[2469]5042       IF (prt_level.ge.9) &
5043            print*,'Attention on met a 0 les thermiques pour phystoke'
[2692]5044       CALL phystokenc ( &
[2469]5045            nlon,klev,pdtphys,longitude_deg,latitude_deg, &
5046            t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
5047            fm_therm,entr_therm, &
5048            cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
5049            frac_impa, frac_nucl, &
[3435]5050            pphis,cell_area,phys_tstep,itap, &
[2469]5051            qx(:,:,ivap),da,phi,mp,upwd,dnwd)
[524]5052
5053
[2469]5054    ENDIF
[524]5055
[2469]5056    !
5057    ! Calculer le transport de l'eau et de l'energie (diagnostique)
5058    !
[4229]5059    CALL transp (paprs,zxtsol, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
5060                 ue, ve, uq, vq, uwat, vwat)
[2469]5061    !
5062    !IM global posePB BEG
5063    IF(1.EQ.0) THEN
5064       !
[4229]5065       CALL transp_lay (paprs,zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, &
[2469]5066            ve_lay, vq_lay, ue_lay, uq_lay)
5067       !
5068    ENDIF !(1.EQ.0) THEN
5069    !IM global posePB END
[4229]5070    !
[2469]5071    ! Accumuler les variables a stocker dans les fichiers histoire:
5072    !
[1279]5073
[2469]5074    !================================================================
5075    ! Conversion of kinetic and potential energy into heat, for
5076    ! parameterisation of subgrid-scale motions
5077    !================================================================
[1753]5078
[2469]5079    d_t_ec(:,:)=0.
5080    forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
[4232]5081    CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx,ivap,iliq,isol, &
[2851]5082         u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
[2469]5083         zmasse,exner,d_t_ec)
5084    t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
[1753]5085
[2469]5086    !=======================================================================
5087    !   SORTIES
5088    !=======================================================================
5089    !
5090    !IM initialisation + calculs divers diag AMIP2
5091    !
5092    include "calcul_divers.h"
5093    !
5094    !IM Interpolation sur les niveaux de pression du NMC
5095    !   -------------------------------------------------
5096    !
5097    include "calcul_STDlev.h"
5098    !
5099    ! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
5100    CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
5101    !
[2496]5102    !cc prw  = eau precipitable
5103    !   prlw = colonne eau liquide
5104    !   prlw = colonne eau solide
[2499]5105    prw(:) = 0.
5106    prlw(:) = 0.
5107    prsw(:) = 0.
5108    DO k = 1, klev
5109       prw(:)  = prw(:)  + q_seri(:,k)*zmasse(:,k)
5110       prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
5111       prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
[2469]5112    ENDDO
5113    !
[4389]5114    IF (ANY(type_trac == ['inca','inco'])) THEN
[655]5115#ifdef INCA
[2469]5116       CALL VTe(VTphysiq)
5117       CALL VTb(VTinca)
[959]5118
[2469]5119       CALL chemhook_end ( &
[3435]5120            phys_tstep, &
[2469]5121            pplay, &
5122            t_seri, &
[3872]5123            tr_seri(:,:,1+nqCO2:nbtr), &
[2469]5124            nbtr, &
5125            paprs, &
5126            q_seri, &
5127            cell_area, &
5128            pphi, &
5129            pphis, &
[2832]5130            zx_rh, &
[4127]5131            aps, bps, ap, bp, lafin)
[959]5132
[2469]5133       CALL VTe(VTinca)
5134       CALL VTb(VTphysiq)
[655]5135#endif
[2692]5136    ENDIF
[655]5137
[4389]5138    IF (type_trac == 'repr') THEN
[4140]5139#ifdef REPROBUS
5140        CALL coord_hyb_rep(paprs, pplay, aps, bps, ap, bp, cell_area)
5141#endif
5142    ENDIF
[1753]5143
[2469]5144    !
5145    ! Convertir les incrementations en tendances
5146    !
5147    IF (prt_level .GE.10) THEN
5148       print *,'Convertir les incrementations en tendances '
5149    ENDIF
5150    !
[2692]5151    IF (mydebug) THEN
5152       CALL writefield_phy('u_seri',u_seri,nbp_lev)
5153       CALL writefield_phy('v_seri',v_seri,nbp_lev)
5154       CALL writefield_phy('t_seri',t_seri,nbp_lev)
5155       CALL writefield_phy('q_seri',q_seri,nbp_lev)
5156    ENDIF
[766]5157
[2469]5158    DO k = 1, klev
5159       DO i = 1, klon
[3435]5160          d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / phys_tstep
5161          d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / phys_tstep
5162          d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / phys_tstep
5163          d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / phys_tstep
5164          d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
[2469]5165          !CR: on ajoute le contenu en glace
[4098]5166          IF (nqo >= 3) THEN
[3435]5167             d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
[2692]5168          ENDIF
[4059]5169          !--ice_sursat: nqo=4, on ajoute rneb
[4098]5170          IF (nqo == 4) THEN
[4059]5171             d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
5172          ENDIF
[2469]5173       ENDDO
5174    ENDDO
5175    !
[4367]5176    ! DC: All iterations are cycled if nqtot==nqo, so no nqtot>nqo condition required
5177    itr = 0
5178    DO iq = 1, nqtot
5179       IF(.NOT.tracers(iq)%isInPhysics) CYCLE
5180       itr = itr+1
5181       DO  k = 1, klev
5182          DO  i = 1, klon
5183             d_qx(i,k,iq) = ( tr_seri(i,k,itr) - qx(i,k,iq) ) / phys_tstep
[2469]5184          ENDDO
5185       ENDDO
[4367]5186    ENDDO
[2469]5187    !
5188    !IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
5189    !IM global posePB      include "write_bilKP_ins.h"
5190    !IM global posePB      include "write_bilKP_ave.h"
5191    !
[1412]5192
[2489]5193    !--OB mass fixer
5194    !--profile is corrected to force mass conservation of water
5195    IF (mass_fixer) THEN
5196    qql2(:)=0.0
[2499]5197    DO k = 1, klev
5198      qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
[2489]5199    ENDDO
5200    DO i = 1, klon
5201      !--compute ratio of what q+ql should be with conservation to what it is
5202      corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
5203      DO k = 1, klev
5204        q_seri(i,k) =q_seri(i,k)*corrqql
5205        ql_seri(i,k)=ql_seri(i,k)*corrqql
5206      ENDDO
5207    ENDDO
5208    ENDIF
5209    !--fin mass fixer
5210
[2469]5211    ! Sauvegarder les valeurs de t et q a la fin de la physique:
5212    !
[2499]5213    u_ancien(:,:)  = u_seri(:,:)
5214    v_ancien(:,:)  = v_seri(:,:)
5215    t_ancien(:,:)  = t_seri(:,:)
5216    q_ancien(:,:)  = q_seri(:,:)
5217    ql_ancien(:,:) = ql_seri(:,:)
5218    qs_ancien(:,:) = qs_seri(:,:)
[4059]5219    rneb_ancien(:,:) = rneb_seri(:,:)
[2499]5220    CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
5221    CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
5222    CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
[2469]5223    ! !! RomP >>>
[4056]5224    IF (nqtot > nqo) tr_ancien(:,:,:) = tr_seri(:,:,:)
[2469]5225    ! !! RomP <<<
5226    !==========================================================================
5227    ! Sorties des tendances pour un point particulier
5228    ! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
5229    ! pour le debug
5230    ! La valeur de igout est attribuee plus haut dans le programme
5231    !==========================================================================
[879]5232
[2692]5233    IF (prt_level.ge.1) THEN
[2469]5234       write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
5235       write(lunout,*) &
5236            'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
5237       write(lunout,*) &
5238            nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
5239            pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
5240            pctsrf(igout,is_sic)
5241       write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
[2692]5242       DO k=1,klev
[2469]5243          write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
5244               d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
5245               d_t_eva(igout,k)
[2692]5246       ENDDO
[2469]5247       write(lunout,*) 'cool,heat'
[2692]5248       DO k=1,klev
[2469]5249          write(lunout,*) cool(igout,k),heat(igout,k)
[2692]5250       ENDDO
[879]5251
[2469]5252       !jyg<     (En attendant de statuer sur le sort de d_t_oli)
5253       !jyg!     write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
5254       !jyg!     do k=1,klev
5255       !jyg!        write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
5256       !jyg!             d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
5257       !jyg!     enddo
5258       write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
[2692]5259       DO k=1,klev
[2469]5260          write(lunout,*) d_t_vdf(igout,k), &
5261               d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
[2692]5262       ENDDO
[2469]5263       !>jyg
[879]5264
[2469]5265       write(lunout,*) 'd_ps ',d_ps(igout)
5266       write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
[2692]5267       DO k=1,klev
[2469]5268          write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
5269               d_qx(igout,k,1),d_qx(igout,k,2)
[2692]5270       ENDDO
5271    ENDIF
[879]5272
[2469]5273    !============================================================
5274    !   Calcul de la temperature potentielle
5275    !============================================================
5276    DO k = 1, klev
5277       DO i = 1, klon
5278          !JYG/IM theta en debut du pas de temps
5279          !JYG/IM       theta(i,k)=t(i,k)*(100000./pplay(i,k))**(RD/RCPD)
5280          !JYG/IM theta en fin de pas de temps de physique
5281          theta(i,k)=t_seri(i,k)*(100000./pplay(i,k))**(RD/RCPD)
5282          ! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
5283          !     MPL 20130625
5284          ! fth_fonctions.F90 et parkind1.F90
5285          ! sinon thetal=theta
5286          !       thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
5287          !    :         ql_seri(i,k))
5288          thetal(i,k)=theta(i,k)
5289       ENDDO
5290    ENDDO
5291    !
[879]5292
[2469]5293    ! 22.03.04 BEG
5294    !=============================================================
5295    !   Ecriture des sorties
5296    !=============================================================
[524]5297#ifdef CPP_IOIPSL
5298
[2469]5299    ! Recupere des varibles calcule dans differents modules
5300    ! pour ecriture dans histxxx.nc
[782]5301
[2469]5302    ! Get some variables from module fonte_neige_mod
5303    CALL fonte_neige_get_vars(pctsrf,  &
[2517]5304         zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
[782]5305
[1507]5306
[2469]5307    !=============================================================
5308    ! Separation entre thermiques et non thermiques dans les sorties
5309    ! de fisrtilp
5310    !=============================================================
[1507]5311
[2692]5312    IF (iflag_thermals>=1) THEN
[2469]5313       d_t_lscth=0.
5314       d_t_lscst=0.
5315       d_q_lscth=0.
5316       d_q_lscst=0.
[2692]5317       DO k=1,klev
5318          DO i=1,klon
5319             IF (ptconvth(i,k)) THEN
[2469]5320                d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5321                d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
[2692]5322             ELSE
[2469]5323                d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5324                d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
[2692]5325             ENDIF
5326          ENDDO
5327       ENDDO
[1507]5328
[2692]5329       DO i=1,klon
[2469]5330          plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
5331          plul_th(i)=prfl(i,1)+psfl(i,1)
[2692]5332       ENDDO
5333    ENDIF
[909]5334
[2469]5335    !On effectue les sorties:
[1791]5336
[2630]5337#ifdef CPP_Dust
5338  CALL phys_output_write_spl(itap, pdtphys, paprs, pphis,  &
5339       pplay, lmax_th, aerosol_couple,                 &
[3776]5340       ok_ade, ok_aie, ivap, ok_sync,                  &
[2630]5341       ptconv, read_climoz, clevSTD,                   &
5342       ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse,      &
5343       flag_aerosol, flag_aerosol_strat, ok_cdnc)
5344#else
[2469]5345    CALL phys_output_write(itap, pdtphys, paprs, pphis,  &
5346         pplay, lmax_th, aerosol_couple,                 &
[3630]5347         ok_ade, ok_aie, ok_volcan, ivap, iliq, isol,    &
[2496]5348         ok_sync, ptconv, read_climoz, clevSTD,          &
[2665]5349         ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
[2469]5350         flag_aerosol, flag_aerosol_strat, ok_cdnc)
[2630]5351#endif
[1791]5352
[2651]5353#ifndef CPP_XIOS
[2590]5354    CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
[2651]5355#endif
[687]5356
[524]5357#endif
5358
[2469]5359    !====================================================================
5360    ! Arret du modele apres hgardfou en cas de detection d'un
5361    ! plantage par hgardfou
5362    !====================================================================
[2235]5363
5364    IF (abortphy==1) THEN
5365       abort_message ='Plantage hgardfou'
[2311]5366       CALL abort_physic (modname,abort_message,1)
[2235]5367    ENDIF
5368
[2469]5369    ! 22.03.04 END
5370    !
5371    !====================================================================
5372    ! Si c'est la fin, il faut conserver l'etat de redemarrage
5373    !====================================================================
5374    !
[782]5375
[3981]5376    ! Disabling calls to the prt_alerte function
5377    alert_first_call = .FALSE.
5378   
[2469]5379    IF (lafin) THEN
5380       itau_phy = itau_phy + itap
5381       CALL phyredem ("restartphy.nc")
5382       !         open(97,form="unformatted",file="finbin")
5383       !         write(97) u_seri,v_seri,t_seri,q_seri
5384       !         close(97)
[3435]5385     
5386       IF (is_omp_master) THEN
5387       
5388         IF (read_climoz >= 1) THEN
5389           IF (is_mpi_root) CALL nf95_close(ncid_climoz)
5390            DEALLOCATE(press_edg_climoz) ! pointer
5391            DEALLOCATE(press_cen_climoz) ! pointer
5392         ENDIF
5393       
[2692]5394       ENDIF
[3435]5395#ifdef CPP_XIOS
5396       IF (is_omp_master) CALL xios_context_finalize
[4127]5397
5398#ifdef INCA
[4389]5399       if (type_trac == 'inca') then
[4127]5400          IF (is_omp_master .and. grid_type==unstructured) THEN
5401             CALL finalize_inca
5402          ENDIF
5403       endif
[3435]5404#endif
[4127]5405
5406#endif
[3461]5407       WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
[2469]5408    ENDIF
[1863]5409
[2469]5410    !      first=.false.
[1863]5411
[2469]5412  END SUBROUTINE physiq
[2418]5413
[2902]5414END MODULE physiq_mod
Note: See TracBrowser for help on using the repository browser.