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

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