source: lmdz_wrf/trunk/WRFV3/lmdz/physiq.F90 @ 1554

Last change on this file since 1554 was 186, checked in by lfita, 10 years ago

Removing checking printings

File size: 165.2 KB
Line 
1! $Id: physiq.F 1795 2013-07-18 08:20:28Z emillour $
2!c#define IO_DEBUG
3
4      SUBROUTINE physiq (nlon,nlev,                                                  &
5       &            debut,lafin,jD_cur, jH_cur,pdtphys,                              &
6       &            paprs,pplay,pphi,pphis,presnivs,clesphy0,                        &
7       &            u,v,t,qx,                                                        &
8       &            flxmass_w,                                                       &
9       &            d_u, d_v, d_t, d_qx, d_ps                                        &
10       &            , dudyn                                                          &
11       &            , PVteta)
12
13      USE ioipsl, only: histbeg, histvert, histdef, histend, histsync,               &
14       &     histwrite, ju2ymds, ymds2ju, ioget_year_len
15      USE comgeomphy
16      USE phys_cal_mod
17      USE write_field_phy
18      USE dimphy
19      USE infotrac
20      USE mod_grid_phy_lmdz
21      USE mod_phys_lmdz_para
22      USE iophy
23      USE misc_mod, mydebug=>debug
24      USE vampir
25!L. Fita, LMD. November 2013
26!      USE pbl_surface_mod, ONLY : pbl_surface
27!      USE pbl_surface_mod, ONLY : pbl_surface, pbl_surface_init
28      USE pbl_surface_mod
29      USE change_srf_frac_mod
30      USE surface_data,     ONLY : type_ocean, ok_veget
31      USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
32      USE phys_state_var_mod ! Variables sauvegardees de la physique
33      USE phys_output_var_mod ! Variables pour les ecritures des sorties
34      USE fonte_neige_mod, ONLY  : fonte_neige_get_vars, fonte_neige_init
35      USE phys_output_mod
36      USE phys_output_ctrlout_mod
37      USE iophy
38      use open_climoz_m, only: open_climoz ! ozone climatology from a file
39      use regr_pr_av_m, only: regr_pr_av
40      use netcdf95, only: nf95_close
41!IM for NMC files
42!c     use netcdf, only: nf90_fill_real
43      use netcdf
44      use mod_phys_lmdz_mpi_data, only: is_mpi_root
45      USE aero_mod
46      use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
47      use conf_phys_m, only: conf_phys
48      use radlwsw_m, only: radlwsw
49      USE control_mod
50#ifdef REPROBUS
51      USE CHEM_REP, ONLY : Init_chem_rep_xjour
52#endif
53      USE indice_sol_mod
54
55!IM stations CFMIP
56      USE CFMIP_point_locations
57
58!L. Fita. July 2013
59      USE orografi
60      USE orografi_strato_mod
61      USE diagphy_mod
62!      USE orbite_mod, ONLY : angle
63      USE phyaqua_mod, ONLY: zenang_an
64      USE nuage_mod, ONLY : diagcld1
65
66! L. Fita, LMD. November 2013.
67!!      USE wrf_lmdz_mod
68!!      USE module_lmdz_variables, ONLY: neige_initialize, limit_initialize
69!!      USE physiq_limit_variables_mod
70!!      USE physiq_limit_variables_mod
71      USE lmdz_wrf_variables_mod
72! L. Fita, LMD. January 2014
73!   Defining with SAVE statement all that variables that should appear in the output
74      USE output_lmdz_NOmodule
75
76      IMPLICIT none
77!>======================================================================
78!!
79!! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
80!!
81!! Objet: Moniteur general de la physique du modele
82!!AA      Modifications quant aux traceurs :
83!!AA                  -  uniformisation des parametrisations ds phytrac
84!!AA                  -  stockage des moyennes des champs necessaires
85!!AA                     en mode traceur off-line
86!!======================================================================
87!!   CLEFS CPP POUR LES IO
88!!   =====================
89#define histNMC
90!c#define histISCCP
91!!======================================================================
92!!    modif   ( P. Le Van ,  12/10/98 )
93!!
94!!  Arguments:
95!!
96!! nlon----input-I-nombre de points horizontaux
97!! nlev----input-I-nombre de couches verticales, doit etre egale a klev
98!! debut---input-L-variable logique indiquant le premier passage
99!! lafin---input-L-variable logique indiquant le dernier passage
100!! jD_cur       -R-jour courant a l'appel de la physique (jour julien)
101!! jH_cur       -R-heure courante a l'appel de la physique (jour julien)
102!! pdtphys-input-R-pas d'integration pour la physique (seconde)
103!! paprs---input-R-pression pour chaque inter-couche (en Pa)
104!! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
105!! pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
106!! pphis---input-R-geopotentiel du sol
107!! presnivs-input_R_pressions approximat. des milieux couches ( en PA)
108!! u-------input-R-vitesse dans la direction X (de O a E) en m/s
109!! v-------input-R-vitesse Y (de S a N) en m/s
110!! t-------input-R-temperature (K)
111!! qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
112!! d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
113!! d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
114!! flxmass_w -input-R- flux de masse verticale
115!! d_u-----output-R-tendance physique de "u" (m/s/s)
116!! d_v-----output-R-tendance physique de "v" (m/s/s)
117!! d_t-----output-R-tendance physique de "t" (K/s)
118!! d_qx----output-R-tendance physique de "qx" (kg/kg/s)
119!! d_ps----output-R-tendance physique de la pression au sol
120!!IM
121!! PVteta--output-R-vorticite potentielle a des thetas constantes
122!!======================================================================
123#include "dimensions.h"
124      integer jjmp1
125      parameter (jjmp1=jjm+1-1/jjm)
126      integer iip1
127      parameter (iip1=iim+1)
128
129#include "regdim.h"
130#include "dimsoil.h"
131#include "clesphys.h"
132#include "temps.h"
133#include "iniprint.h"
134#include "thermcell.h"
135!c======================================================================
136      LOGICAL ok_cvl  ! pour activer le nouveau driver pour convection KE
137      PARAMETER (ok_cvl=.TRUE.)
138      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
139      PARAMETER (ok_gust=.FALSE.)
140      integer iflag_radia     ! active ou non le rayonnement (MPL)
141      save iflag_radia
142!$OMP THREADPRIVATE(iflag_radia)
143!c======================================================================
144      LOGICAL check ! Verifier la conservation du modele en eau
145      PARAMETER (check=.FALSE.)
146      LOGICAL ok_stratus ! Ajouter artificiellement les stratus
147      PARAMETER (ok_stratus=.FALSE.)
148!c======================================================================
149      REAL amn, amx
150      INTEGER igout
151!c======================================================================
152!c Clef controlant l'activation du cycle diurne:
153!ccc      LOGICAL cycle_diurne
154!ccc      PARAMETER (cycle_diurne=.FALSE.)
155!c======================================================================
156!c Modele thermique du sol, a activer pour le cycle diurne:
157!ccc      LOGICAL soil_model
158!ccc      PARAMETER (soil_model=.FALSE.)
159!c======================================================================
160!c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
161!c le calcul du rayonnement est celle apres la precipitation des nuages.
162!c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
163!c la condensation et la precipitation. Cette cle augmente les impacts
164!c radiatifs des nuages.
165!ccc      LOGICAL new_oliq
166!ccc      PARAMETER (new_oliq=.FALSE.)
167!c======================================================================
168!c Clefs controlant deux parametrisations de l'orographie:
169!cc      LOGICAL ok_orodr
170!ccc      PARAMETER (ok_orodr=.FALSE.)
171!ccc      LOGICAL ok_orolf
172!ccc      PARAMETER (ok_orolf=.FALSE.)
173!c======================================================================
174      LOGICAL ok_journe ! sortir le fichier journalier
175      save ok_journe
176!$OMP THREADPRIVATE(ok_journe)
177!c
178      LOGICAL ok_mensuel ! sortir le fichier mensuel
179      save ok_mensuel
180!$OMP THREADPRIVATE(ok_mensuel)
181!c
182      LOGICAL ok_instan ! sortir le fichier instantane
183      save ok_instan
184!$OMP THREADPRIVATE(ok_instan)
185!c
186      LOGICAL ok_LES ! sortir le fichier LES
187      save ok_LES                           
188!$OMP THREADPRIVATE(ok_LES)                 
189!c
190      LOGICAL callstats ! sortir le fichier stats
191      save callstats                           
192!$OMP THREADPRIVATE(callstats)                 
193!c
194      LOGICAL ok_region ! sortir le fichier regional
195      PARAMETER (ok_region=.FALSE.)
196!c======================================================================
197! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
198!      real weak_inversion(klon),dthmin(klon)
199      real seuil_inversion
200      save seuil_inversion
201!$OMP THREADPRIVATE(seuil_inversion)
202      integer iflag_ratqs
203      save iflag_ratqs
204!$OMP THREADPRIVATE(iflag_ratqs)
205      real facteur
206
207      REAL zz,znum,zden
208      REAL wmax_th(klon)
209! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
210!      REAL zmax_th(klon)
211      REAL tau_overturning_th(klon)
212
213! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
214!      integer lmax_th(klon)
215      integer limbas(klon)
216      real ratqscth(klon,klev)
217      real ratqsdiff(klon,klev)
218      real zqsatth(klon,klev)
219
220!c======================================================================
221!c
222      INTEGER ivap          ! indice de traceurs pour vapeur d'eau
223      PARAMETER (ivap=1)
224      INTEGER iliq          ! indice de traceurs pour eau liquide
225      PARAMETER (iliq=2)
226
227!c
228!c
229!c Variables argument:
230!c
231      INTEGER nlon
232      INTEGER nlev
233      REAL, intent(in):: jD_cur, jH_cur
234
235      REAL pdtphys
236      LOGICAL debut, lafin
237      REAL paprs(klon,klev+1)
238      REAL pplay(klon,klev)
239      REAL pphi(klon,klev)
240      REAL pphis(klon)
241      REAL presnivs(klev)
242      REAL znivsig(klev)
243      real pir
244
245      REAL u(klon,klev)
246      REAL v(klon,klev)
247! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
248!      REAL t(klon,klev),theta(klon,klev),thetal(klon,klev)
249      REAL t(klon,klev),thetal(klon,klev)
250
251!c thetal: ligne suivante a decommenter si vous avez les fichiers     MPL 20130625
252!c fth_fonctions.F90 et parkind1.F90
253!c sinon thetal=theta
254!c     REAL fth_thetae,fth_thetav,fth_thetal
255      REAL qx(klon,klev,nqtot)
256      REAL flxmass_w(klon,klev)
257! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
258!      REAL omega(klon,klev) ! vitesse verticale en Pa/s
259      REAL d_u(klon,klev)
260      REAL d_v(klon,klev)
261      REAL d_t(klon,klev)
262      REAL d_qx(klon,klev,nqtot)
263      REAL d_ps(klon)
264! Variables pour le transport convectif
265      real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
266! Variables pour le lessivage convectif
267! RomP >>>
268      real phi2(klon,klev,klev)
269      real d1a(klon,klev),dam(klon,klev)
270      real ev(klon,klev),ep(klon,klev)
271      real clw(klon,klev),elij(klon,klev,klev)
272      real epmlmMm(klon,klev,klev),eplaMm(klon,klev)
273! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
274!      real wdtrainA(klon,klev),wdtrainM(klon,klev)
275
276! RomP <<<
277!IM definition dynamique o_trac dans phys_output_open
278!      type(ctrl_out) :: o_trac(nqtot)
279!c
280!IM Amip2 PV a theta constante
281!c
282      INTEGER nbteta
283      PARAMETER(nbteta=3)
284      CHARACTER*3 ctetaSTD(nbteta)
285      DATA ctetaSTD/'350','380','405'/
286      SAVE ctetaSTD
287!$OMP THREADPRIVATE(ctetaSTD)
288      REAL rtetaSTD(nbteta)
289      DATA rtetaSTD/350., 380., 405./
290      SAVE rtetaSTD
291!$OMP THREADPRIVATE(rtetaSTD)     
292!c
293      REAL PVteta(klon,nbteta)
294      REAL zx_tmp_3dte(iim,jjmp1,nbteta)
295!c
296!MI Amip2 PV a theta constante
297
298!cym      INTEGER klevp1, klevm1
299!cym      PARAMETER(klevp1=klev+1,klevm1=klev-1)
300!cym#include "raddim.h"
301!c
302!c
303!IM Amip2
304!c variables a une pression donnee
305!c
306#include "declare_STDlev.h"
307!c
308      CHARACTER*4 bb2
309      CHARACTER*2 bb3
310!c
311#include "radopt.h"
312!c
313!c
314!c prw: precipitable water
315! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
316!      real prw(klon)
317
318      REAL convliq(klon,klev)  ! eau liquide nuageuse convective
319      REAL convfra(klon,klev)  ! fraction nuageuse convective
320
321      REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut
322      REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree
323      REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut
324      REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
325
326      INTEGER linv, kp1
327!c flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
328!c flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
329! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
330!      REAL flwp(klon), fiwp(klon)
331!      REAL flwc(klon,klev), fiwc(klon,klev)
332      REAL flwp_c(klon), fiwp_c(klon)
333      REAL flwc_c(klon,klev), fiwc_c(klon,klev)
334      REAL flwp_s(klon), fiwp_s(klon)
335      REAL flwc_s(klon,klev), fiwc_s(klon,klev)
336
337! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
338!      REAL evap_pot(klon,nbsrf)
339
340!IM ISCCP simulator v3.4
341!c dans clesphys.h top_height, overlap
342!cv3.4
343      INTEGER debug, debugcol
344!cym      INTEGER npoints
345!cym      PARAMETER(npoints=klon)
346!c
347      INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night
348      INTEGER nregISCtot
349      PARAMETER(nregISCtot=1)
350!c
351!c imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
352!c y compris pour 1 point
353!c imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
354!c jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
355      INTEGER imin_debut, nbpti
356      INTEGER jmin_debut, nbptj
357!IM parametres ISCCP BEG
358      INTEGER nbapp_isccp
359!     INTEGER nbapp_isccp,isccppas
360!     PARAMETER(isccppas=6) !appel du simulateurs tous les 6pas de temps de la physique
361!                           !i.e. toutes les 3 heures
362      INTEGER n
363      INTEGER ifreq_isccp(napisccp), freqin_pdt(napisccp)
364      DATA ifreq_isccp/3/
365      SAVE ifreq_isccp
366!$OMP THREADPRIVATE(ifreq_isccp)
367      CHARACTER*5 typinout(napisccp)
368      DATA typinout/'i3od'/
369      SAVE typinout
370!$OMP THREADPRIVATE(typinout)
371!IM verif boxptop BEG
372      CHARACTER*1 verticaxe(napisccp)
373      DATA verticaxe/'1'/
374      SAVE verticaxe
375!$OMP THREADPRIVATE(verticaxe)
376!IM verif boxptop END
377      INTEGER nvlev(napisccp)
378!c     INTEGER nvlev
379      REAL t1, aa
380      REAL seed_re(klon,napisccp)
381!cym !!!! A voir plus tard
382!cym      INTEGER iphy(iim,jjmp1)
383!IM parametres ISCCP END
384!c
385!c ncol = nb. de sous-colonnes pour chaque maille du GCM
386!c ncolmx = No. max. de sous-colonnes pour chaque maille du GCM
387!c      INTEGER ncol(napisccp), ncolmx, seed(klon,napisccp)
388      INTEGER,SAVE :: ncol(napisccp)
389!$OMP THREADPRIVATE(ncol)
390      INTEGER ncolmx, seed(klon,napisccp)
391      REAL nbsunlit(nregISCtot,klon,napisccp)  !nbsunlit : moyenne de sunlit
392!c     PARAMETER(ncolmx=1500)
393      PARAMETER(ncolmx=300)
394!c
395!IM verif boxptop BEG
396      REAL vertlev(ncolmx,napisccp)
397!IM verif boxptop END
398!c
399      REAL,SAVE :: tautab_omp(0:255),tautab(0:255)
400      INTEGER,SAVE :: invtau_omp(-20:45000),invtau(-20:45000)
401!$OMP THREADPRIVATE(tautab,invtau)
402      REAL emsfc_lw
403      PARAMETER(emsfc_lw=0.99)
404!c     REAL    ran0                      ! type for random number fuction
405!c
406      REAL cldtot(klon,klev)
407!c variables de haut en bas pour le simulateur ISCCP
408      REAL dtau_s(klon,klev) !tau nuages startiformes
409      REAL dtau_c(klon,klev) !tau nuages convectifs
410      REAL dem_s(klon,klev)  !emissivite nuages startiformes
411      REAL dem_c(klon,klev)  !emissivite nuages convectifs
412!c
413!c variables de haut en bas pour le simulateur ISCCP
414      REAL pfull(klon,klev)
415      REAL phalf(klon,klev+1)
416      REAL qv(klon,klev)
417      REAL cc(klon,klev)
418      REAL conv(klon,klev)
419      REAL dtau_sH2B(klon,klev)
420      REAL dtau_cH2B(klon,klev)
421      REAL at(klon,klev)
422      REAL dem_sH2B(klon,klev)
423      REAL dem_cH2B(klon,klev)
424!c
425      INTEGER kmax, lmax, lmax3
426      PARAMETER(kmax=8, lmax=8, lmax3=3)
427      INTEGER kmaxm1, lmaxm1
428      PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
429      INTEGER iimx7, jjmx7, jjmp1x7
430      PARAMETER(iimx7=iim*kmaxm1, jjmx7=jjm*lmaxm1,                                  &
431       &jjmp1x7=jjmp1*lmaxm1)
432!c
433!c output from ISCCP simulator
434      REAL fq_isccp(klon,kmaxm1,lmaxm1,napisccp)
435      REAL fq_is_true(klon,kmaxm1,lmaxm1,napisccp)
436      REAL totalcldarea(klon,napisccp)
437      REAL meanptop(klon,napisccp)
438      REAL meantaucld(klon,napisccp)
439      REAL boxtau(klon,ncolmx,napisccp)
440      REAL boxptop(klon,ncolmx,napisccp)
441      REAL zx_tmp_fi3d_bx(klon,ncolmx)
442      REAL zx_tmp_3d_bx(iim,jjmp1,ncolmx)
443!c
444      REAL cld_fi3d(klon,lmax3)
445      REAL cld_3d(iim,jjmp1,lmax3)
446!c
447      INTEGER iw, iwmax
448      REAL wmin, pas_w
449!c     PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
450!IM 051005     PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
451      PARAMETER(wmin=-100.,pas_w=10.,iwmax=20)
452      REAL o500(klon)
453!c
454
455!c sorties ISCCP
456
457      integer nid_isccp
458      save nid_isccp       
459!$OMP THREADPRIVATE(nid_isccp)
460
461      REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
462      DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
463      SAVE zx_tau
464      DATA zx_pc/180., 310., 440., 560., 680., 800., 1000./
465      SAVE zx_pc
466!$OMP THREADPRIVATE(zx_tau,zx_pc)
467!c cldtopres pression au sommet des nuages
468      REAL cldtopres(lmaxm1), cldtopres3(lmax3)
469      DATA cldtopres/180., 310., 440., 560., 680., 800., 1000./
470      DATA cldtopres3/440., 680., 1000./
471      SAVE cldtopres,cldtopres3
472!$OMP THREADPRIVATE(cldtopres,cldtopres3)
473!IM 051005 BEG
474      INTEGER komega, nhoriRD
475
476!c taulev: numero du niveau de tau dans les sorties ISCCP
477      CHARACTER *4 taulev(kmaxm1)
478!c     DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
479      DATA taulev/'tau0','tau1','tau2','tau3','tau4','tau5','tau6'/
480      CHARACTER *3 pclev(lmaxm1)
481      DATA pclev/'pc1','pc2','pc3','pc4','pc5','pc6','pc7'/
482      SAVE taulev,pclev
483!$OMP THREADPRIVATE(taulev,pclev)
484!c
485!c cnameisccp
486      CHARACTER *29 cnameisccp(lmaxm1,kmaxm1)
487!IM bad 151205     DATA cnameisccp/'pc< 50hPa, tau< 0.3',
488      DATA cnameisccp/'pc= 50-180hPa, tau< 0.3',                                     &
489       &                'pc= 180-310hPa, tau< 0.3',                                  &
490       &                'pc= 310-440hPa, tau< 0.3',                                  &
491       &                'pc= 440-560hPa, tau< 0.3',                                  &
492       &                'pc= 560-680hPa, tau< 0.3',                                  &
493       &                'pc= 680-800hPa, tau< 0.3',                                  &
494       &                'pc= 800-1000hPa, tau< 0.3',                                 &
495       &                'pc= 50-180hPa, tau= 0.3-1.3',                               &
496       &                'pc= 180-310hPa, tau= 0.3-1.3',                              &
497       &                'pc= 310-440hPa, tau= 0.3-1.3',                              &
498       &                'pc= 440-560hPa, tau= 0.3-1.3',                              &
499       &                'pc= 560-680hPa, tau= 0.3-1.3',                              &
500       &                'pc= 680-800hPa, tau= 0.3-1.3',                              &
501       &                'pc= 800-1000hPa, tau= 0.3-1.3',                             &
502       &                'pc= 50-180hPa, tau= 1.3-3.6',                               &
503       &                'pc= 180-310hPa, tau= 1.3-3.6',                              &
504       &                'pc= 310-440hPa, tau= 1.3-3.6',                              &
505       &                'pc= 440-560hPa, tau= 1.3-3.6',                              &
506       &                'pc= 560-680hPa, tau= 1.3-3.6',                              &
507       &                'pc= 680-800hPa, tau= 1.3-3.6',                              &
508       &                'pc= 800-1000hPa, tau= 1.3-3.6',                             &
509       &                'pc= 50-180hPa, tau= 3.6-9.4',                               &
510       &                'pc= 180-310hPa, tau= 3.6-9.4',                              &
511       &                'pc= 310-440hPa, tau= 3.6-9.4',                              &
512       &                'pc= 440-560hPa, tau= 3.6-9.4',                              &
513       &                'pc= 560-680hPa, tau= 3.6-9.4',                              &
514       &                'pc= 680-800hPa, tau= 3.6-9.4',                              &
515       &                'pc= 800-1000hPa, tau= 3.6-9.4',                             &
516       &                'pc= 50-180hPa, tau= 9.4-23',                                &
517       &                'pc= 180-310hPa, tau= 9.4-23',                               &
518       &                'pc= 310-440hPa, tau= 9.4-23',                               &
519       &                'pc= 440-560hPa, tau= 9.4-23',                               &
520       &                'pc= 560-680hPa, tau= 9.4-23',                               &
521       &                'pc= 680-800hPa, tau= 9.4-23',                               &
522       &                'pc= 800-1000hPa, tau= 9.4-23',                              &
523       &                'pc= 50-180hPa, tau= 23-60',                                 &
524       &                'pc= 180-310hPa, tau= 23-60',                                &
525       &                'pc= 310-440hPa, tau= 23-60',                                &
526       &                'pc= 440-560hPa, tau= 23-60',                                &
527       &                'pc= 560-680hPa, tau= 23-60',                                &
528       &                'pc= 680-800hPa, tau= 23-60',                                &
529       &                'pc= 800-1000hPa, tau= 23-60',                               &
530       &                'pc= 50-180hPa, tau> 60.',                                   &
531       &                'pc= 180-310hPa, tau> 60.',                                  &
532       &                'pc= 310-440hPa, tau> 60.',                                  &
533       &                'pc= 440-560hPa, tau> 60.',                                  &
534       &                'pc= 560-680hPa, tau> 60.',                                  &
535       &                'pc= 680-800hPa, tau> 60.',                                  &
536       &                'pc= 800-1000hPa, tau> 60.'/
537       SAVE cnameisccp
538!$OMP THREADPRIVATE(cnameisccp)
539!c
540!c     REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
541!c     INTEGER nhorix7
542!IM: region='3d' <==> sorties en global
543      CHARACTER*3 region
544      PARAMETER(region='3d')
545!c
546!IM ISCCP simulator v3.4
547!c
548      logical ok_hf
549!c
550      integer nid_hf, nid_hf3d
551      save ok_hf, nid_hf, nid_hf3d
552!$OMP THREADPRIVATE(ok_hf, nid_hf, nid_hf3d)
553!c  QUESTION : noms de variables ?
554
555      INTEGER        longcles
556      PARAMETER    ( longcles = 20 )
557      REAL clesphy0( longcles      )
558!c
559!c Variables propres a la physique
560      INTEGER itap
561      SAVE itap                   ! compteur pour la physique
562!$OMP THREADPRIVATE(itap)
563!c
564! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
565!      real slp(klon) ! sea level pressure
566!c
567! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
568!      REAL fevap(klon,nbsrf)
569!      REAL fluxlat(klon,nbsrf)
570!c
571! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
572!      REAL qsol(klon)
573      REAL,save ::  solarlong0
574!$OMP THREADPRIVATE(solarlong0)
575
576!c
577!c  Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
578!c
579!IM 141004     REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
580      REAL zulow(klon),zvlow(klon)
581!c
582      INTEGER igwd,idx(klon),itest(klon)
583!c
584! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
585!      REAL agesno(klon,nbsrf)
586
587!c
588!c      REAL,allocatable,save :: run_off_lic_0(:)
589!cc$OMP THREADPRIVATE(run_off_lic_0)
590!cym      SAVE run_off_lic_0
591!cKE43
592!c Variables liees a la convection de K. Emanuel (sb):
593!c
594      REAL bas, top             ! cloud base and top levels
595      SAVE bas
596      SAVE top
597!$OMP THREADPRIVATE(bas, top)
598
599      REAL wdn(klon), tdn(klon), qdn(klon)
600!c
601!c=================================================================================================
602!cCR04.12.07: on ajoute les nouvelles variables du nouveau schema de convection avec poches froides
603!c Variables li\'ees \`a la poche froide (jyg)
604
605      REAL mip(klon,klev)  ! mass flux shed by the adiab ascent at each level
606! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
607!      REAL Vprecip(klon,klev+1)   ! precipitation vertical profile
608!c
609      REAL wape_prescr, fip_prescr
610      INTEGER it_wape_prescr
611      SAVE wape_prescr, fip_prescr, it_wape_prescr
612!$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
613!c
614!c variables supplementaires de concvl
615      REAL Tconv(klon,klev)
616      REAL ment(klon,klev,klev),sij(klon,klev,klev)
617      REAL dd_t(klon,klev),dd_q(klon,klev)
618
619      real, save :: alp_bl_prescr=0.
620      real, save :: ale_bl_prescr=0.
621
622      real, save :: ale_max=1000.
623      real, save :: alp_max=2.
624
625      real, save :: wake_s_min_lsp=0.1
626
627!$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
628!$OMP THREADPRIVATE(ale_max,alp_max)
629!$OMP THREADPRIVATE(wake_s_min_lsp)
630
631! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
632!      real ale_wake(klon)
633!      real alp_wake(klon)
634
635      real ok_wk_lsp(klon)
636
637!cRC
638!c Variables li\'ees \`a la poche froide (jyg et rr)
639!c Version diagnostique pour l'instant : pas de r\'etroaction sur la convection
640
641      REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
642
643      REAL wake_dth(klon,klev)        ! wake : temp pot difference
644
645      REAL wake_d_deltat_gw(klon,klev)! wake : delta T tendency due to Gravity Wave (/s)
646      REAL wake_omgbdth(klon,klev)    ! Wake : flux of Delta_Theta transported by LS omega
647      REAL wake_dp_omgb(klon,klev)    ! Wake : vertical gradient of large scale omega
648      REAL wake_dtKE(klon,klev)       ! Wake : differential heating (wake - unpertubed) CONV
649      REAL wake_dqKE(klon,klev)       ! Wake : differential moistening (wake - unpertubed) CONV
650      REAL wake_dtPBL(klon,klev)      ! Wake : differential heating (wake - unpertubed) PBL
651      REAL wake_dqPBL(klon,klev)      ! Wake : differential moistening (wake - unpertubed) PBL
652! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
653!      REAL wake_omg(klon,klev)        ! Wake : velocity difference (wake - unpertubed)
654      REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
655      REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
656      REAL wake_spread(klon,klev)     ! spreading term in wake_delt
657!c
658!cpourquoi y'a pas de save??
659! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
660!      REAL wake_h(klon)               ! Wake : hauteur de la poche froide
661!c
662      INTEGER wake_k(klon)            ! Wake sommet
663!c
664      REAL t_undi(klon,klev)               ! temperature moyenne dans la zone non perturbee
665      REAL q_undi(klon,klev)               ! humidite moyenne dans la zone non perturbee
666!c
667!cjyg
668!ccc      REAL wake_pe(klon)              ! Wake potential energy - WAPE
669
670      REAL wake_gfl(klon)             ! Gust Front Length
671      REAL wake_dens(klon)
672!c
673!c
674      REAL dt_dwn(klon,klev)
675      REAL dq_dwn(klon,klev)
676      REAL wdt_PBL(klon,klev)
677      REAL udt_PBL(klon,klev)
678      REAL wdq_PBL(klon,klev)
679      REAL udq_PBL(klon,klev)
680      REAL M_dwn(klon,klev)
681      REAL M_up(klon,klev)
682      REAL dt_a(klon,klev)
683      REAL dq_a(klon,klev)
684      REAL, SAVE :: alp_offset
685!$OMP THREADPRIVATE(alp_offset)
686
687!c
688!cRR:fin declarations poches froides
689!c=======================================================================================================
690
691      REAL zw2(klon,klev+1)
692
693! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
694!      REAL fraca(klon,klev+1)       
695      REAL ztv(klon,klev),ztva(klon,klev)
696      REAL zpspsk(klon,klev)
697      REAL ztla(klon,klev),zqla(klon,klev)
698      REAL zthl(klon,klev)
699
700!ccc nrlmd le 10/04/2012
701
702!c--------Stochastic Boundary Layer Triggering: ALE_BL--------
703!c---Propri\'et\'es du thermiques au LCL
704      real zlcl_th(klon)                                     ! Altitude du LCL calcul\'e continument (pcon dans thermcell_main.F90)
705      real fraca0(klon)                                      ! Fraction des thermiques au LCL
706      real w0(klon)                                          ! Vitesse des thermiques au LCL
707      real w_conv(klon)                                      ! Vitesse verticale de grande \'echelle au LCL
708      real tke0(klon,klev+1)                                 ! TKE au début du pas de temps
709      real therm_tke_max0(klon)                              ! TKE dans les thermiques au LCL
710      real env_tke_max0(klon)                                ! TKE dans l'environnement au LCL
711
712!c---Spectre de thermiques de type 2 au LCL
713      real n2(klon),s2(klon)
714      real ale_bl_stat(klon)
715
716!c---D\'eclenchement stochastique
717      integer :: tau_trig(klon)
718      real proba_notrig(klon)
719      real random_notrig(klon)
720
721!c--------Statistical Boundary Layer Closure: ALP_BL--------
722!c---Profils de TKE dans et hors du thermique
723      real pbl_tke_input(klon,klev+1,nbsrf)
724      real therm_tke_max(klon,klev)                          ! Profil de TKE dans les thermiques
725      real env_tke_max(klon,klev)                            ! Profil de TKE dans l'environnement
726
727!c---Fermeture statistique
728      real alp_bl_det(klon)                                     ! ALP d\'terministe du thermique unique
729      real alp_bl_fluct_m(klon)                                 ! ALP li\'ee aux fluctuations de flux de masse sous-nuageux
730      real alp_bl_fluct_tke(klon)                               ! ALP li\'ee aux fluctuations d'\'energie cin\'etique sous-nuageuse
731      real alp_bl_conv(klon)                                    ! ALP li\'ee \`a grande \'echelle
732      real alp_bl_stat(klon)                                    ! ALP totale
733
734!ccc fin nrlmd le 10/04/2012
735
736!c Variables locales pour la couche limite (al1):
737!c
738!cAl1      REAL pblh(klon)           ! Hauteur de couche limite
739!cAl1      SAVE pblh
740!c34EK
741!c
742!c Variables locales:
743!c
744! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
745!      REAL cdragh(klon) ! drag coefficient pour T and Q
746!      REAL cdragm(klon) ! drag coefficient pour vent
747!cAA
748!cAA  Pour phytrac
749      REAL u1(klon)             ! vents dans la premiere couche U
750      REAL v1(klon)             ! vents dans la premiere couche V
751
752! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule\n
753!      REAL zxffonte(klon), zxfqcalving(klon),zxfqfonte(klon)
754
755!@$$      LOGICAL offline           ! Controle du stockage ds "physique"
756!@$$      PARAMETER (offline=.false.)
757!@$$      INTEGER physid
758      REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
759      REAL frac_nucl(klon,klev) ! idem (nucleation)
760! RomP >>>
761      REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
762! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
763!      REAL beta_prec(klon,klev)       ! taux de conv de l'eau cond (utilise)
764! RomP <<<
765      INTEGER       :: iii
766      REAL          :: calday
767
768!IM cf FH pour Tiedtke 080604
769      REAL rain_tiedtke(klon),snow_tiedtke(klon)
770!c
771!IM 050204 END
772! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
773!      REAL evap(klon), devap(klon) ! evaporation et sa derivee
774      REAL devap(klon) ! evaporation et sa derivee
775
776! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
777!      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
778      REAL dsens(klon) ! chaleur sensible et sa derivee
779
780! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
781!      REAL bils(klon) ! bilan de chaleur au sol
782!      REAL wfbilo(klon,nbsrf) ! bilan d'eau, pour chaque
783!C                             ! type de sous-surface et pondere par la fraction
784!      REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
785!C                             ! type de sous-surface et pondere par la fraction
786      REAL slab_wfbils(klon)  ! bilan de chaleur au sol pour le cas de slab, sur les points d'ocean
787
788! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
789!      REAL fder(klon)         
790!      REAL ve(klon) ! integr. verticale du transport meri. de l'energie
791!      REAL vq(klon) ! integr. verticale du transport meri. de l'eau
792!      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
793!      REAL uq(klon) ! integr. verticale du transport zonal de l'eau
794!c
795
796! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
797!      REAL frugs(klon,nbsrf)
798!      REAL zxrugs(klon) ! longueur de rugosite
799!c
800!c Conditions aux limites
801!c
802!
803      REAL :: day_since_equinox
804! Date de l'equinoxe de printemps
805      INTEGER, parameter :: mth_eq=3, day_eq=21
806      REAL :: jD_eq
807
808      LOGICAL, parameter :: new_orbit = .true.
809
810!c
811      INTEGER lmt_pas
812      SAVE lmt_pas                ! frequence de mise a jour
813!$OMP THREADPRIVATE(lmt_pas)
814! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
815!      real zmasse(klon, llm),exner(klon, llm)
816      real exner(klon, llm)
817!C     (column-density of mass of air in a cell, in kg m-2)
818      real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
819
820!IM sorties
821      REAL un_jour
822      PARAMETER(un_jour=86400.)
823!c======================================================================
824!c
825!c Declaration des procedures appelees
826!c
827      EXTERNAL angle     ! calculer angle zenithal du soleil
828      EXTERNAL alboc     ! calculer l'albedo sur ocean
829      EXTERNAL ajsec     ! ajustement sec
830      EXTERNAL conlmd    ! convection (schema LMD)
831!cKE43
832      EXTERNAL conema3  ! convect4.3
833      EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)
834!cAA
835      EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
836!c                          ! stockage des coefficients necessaires au
837!c                          ! lessivage OFF-LINE et ON-LINE
838      EXTERNAL hgardfou  ! verifier les temperatures
839      EXTERNAL nuage     ! calculer les proprietes radiatives
840!CC      EXTERNAL o3cm      ! initialiser l'ozone
841      EXTERNAL orbite    ! calculer l'orbite terrestre
842      EXTERNAL phyetat0  ! lire l'etat initial de la physique
843      EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
844      EXTERNAL suphel    ! initialiser certaines constantes
845      EXTERNAL transp    ! transport total de l'eau et de l'energie
846      EXTERNAL ecribina  ! ecrire le fichier binaire global
847      EXTERNAL ecribins  ! ecrire le fichier binaire global
848      EXTERNAL ecrirega  ! ecrire le fichier binaire regional
849      EXTERNAL ecriregs  ! ecrire le fichier binaire regional
850!IM
851      EXTERNAL haut2bas  !variables de haut en bas
852      EXTERNAL ini_undefSTD  !initialise a 0 une variable a 1 niveau de pression
853      EXTERNAL undefSTD      !somme les valeurs definies d'1 var a 1 niveau de pression
854!c     EXTERNAL moy_undefSTD  !moyenne d'1 var a 1 niveau de pression
855!c     EXTERNAL moyglo_aire   !moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
856!c                            !par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
857!c
858!c Variables locales
859!c
860      REAL rhcl(klon,klev)    ! humiditi relative ciel clair
861      REAL dialiq(klon,klev)  ! eau liquide nuageuse
862      REAL diafra(klon,klev)  ! fraction nuageuse
863      REAL cldliq(klon,klev)  ! eau liquide nuageuse
864! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
865!      REAL cldfra(klon,klev)  ! fraction nuageuse
866!      REAL cldtau(klon,klev)  ! epaisseur optique
867!      REAL cldemi(klon,klev)  ! emissivite infrarouge
868!c
869!CXXX PB
870! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
871      REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
872!      REAL fluxt(klon,klev, nbsrf)   ! flux turbulent de chaleur
873!      REAL fluxu(klon,klev, nbsrf)   ! flux turbulent de vitesse u
874!      REAL fluxv(klon,klev, nbsrf)   ! flux turbulent de vitesse v
875!c
876      REAL zxfluxt(klon, klev)
877      REAL zxfluxq(klon, klev)
878      REAL zxfluxu(klon, klev)
879      REAL zxfluxv(klon, klev)
880!CXXX
881!c
882! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
883!      REAL fsollw(klon, nbsrf)   ! bilan flux IR pour chaque sous surface
884!      REAL fsolsw(klon, nbsrf)   ! flux solaire absorb. pour chaque sous surface
885!c Le rayonnement n'est pas calcule tous les pas, il faut donc
886!c                      sauvegarder les sorties du rayonnement
887!cym      SAVE  heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
888!cym      SAVE  sollwdownclr, toplwdown, toplwdownclr
889!cym      SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
890!c
891      INTEGER itaprad
892      SAVE itaprad
893!$OMP THREADPRIVATE(itaprad)
894!c
895      REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
896      REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
897!c
898! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
899!      REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
900!      REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
901!c
902!      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
903      REAL zxsnow_dummy(klon)
904!c
905      REAL dist, rmu0(klon), fract(klon)
906      REAL zdtime, zlongi
907!c
908      CHARACTER*2 str2
909      CHARACTER*2 iqn
910!c
911      REAL qcheck
912      REAL z_avant(klon), z_apres(klon), z_factor(klon)
913      LOGICAL zx_ajustq
914!c
915      REAL za, zb
916      REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
917      real zqsat(klon,klev)
918      INTEGER i, k, iq, ig, j, nsrf, ll, l, iiq, iff
919      REAL t_coup
920      PARAMETER (t_coup=234.0)
921!c
922! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
923!      REAL zphi(klon,klev)
924!cym A voir plus tard !!
925!cym      REAL zx_relief(iim,jjmp1)
926!cym      REAL zx_aire(iim,jjmp1)
927!c
928!c Grandeurs de sorties
929! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
930!      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
931!      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
932!      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
933      REAL s_capCL(klon)
934      REAL s_oliqCL(klon), s_cteiCL(klon)
935      REAL s_trmb1(klon), s_trmb2(klon)
936      REAL s_trmb3(klon)
937!cKE43
938!c Variables locales pour la convection de K. Emanuel (sb):
939!c
940! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
941!      REAL upwd(klon,klev)      ! saturated updraft mass flux
942!      REAL dnwd(klon,klev)      ! saturated downdraft mass flux
943!      REAL dnwd0(klon,klev)     ! unsaturated downdraft mass flux
944      REAL tvp(klon,klev)       ! virtual temp of lifted parcel
945!      REAL plcl(klon)           ! Lifting Condensation Level
946!      REAL plfc(klon)           ! Level of Free Convection
947!      REAL wbeff(klon)          ! saturated updraft velocity at LFC
948      CHARACTER*40 capemaxcels  !max(CAPE)
949
950      REAL rflag(klon)          ! flag fonctionnement de convect
951      INTEGER iflagctrl(klon)          ! flag fonctionnement de convect
952
953!c -- convect43:
954      INTEGER ntra              ! nb traceurs pour convect4.3
955      REAL pori_con(klon)    ! pressure at the origin level of lifted parcel
956      REAL dtma_con(klon),dtlcl_con(klon)
957      REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
958      REAL dplcldt(klon), dplcldr(klon)
959!c?     .     condm_con(klon,klev),conda_con(klon,klev),
960!c?     .     mr_con(klon,klev),ep_con(klon,klev)
961!c?     .    ,sadiab(klon,klev),wadiab(klon,klev)
962!c --
963!c34EK
964!c
965!c Variables du changement
966!c
967!c con: convection
968!c lsc: condensation a grande echelle (Large-Scale-Condensation)
969!c ajs: ajustement sec
970!c eva: evaporation de l'eau liquide nuageuse
971!c vdf: couche limite (Vertical DiFfusion)
972! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
973!      REAL rneb(klon,klev)
974
975! tendance nulles
976      REAL du0(klon,klev),dv0(klon,klev),dq0(klon,klev),dql0(klon,klev)
977
978!c
979!*********************************************************
980!*     declarations
981     
982!*********************************************************
983!IM 081204 END
984!c
985! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
986!      REAL pmfu(klon,klev), pmfd(klon,klev)
987      REAL pen_u(klon,klev), pen_d(klon,klev)
988      REAL pde_u(klon,klev), pde_d(klon,klev)
989      INTEGER kcbot(klon), kctop(klon), kdtop(klon)
990! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
991!      REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
992!      REAL prfl(klon,klev+1), psfl(klon,klev+1)
993!c
994! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
995!      REAL rain_lsc(klon)
996!      REAL snow_lsc(klon)
997!c
998      REAL ratqsc(klon,klev)
999      real ratqsbas,ratqshaut,tau_ratqs
1000      save ratqsbas,ratqshaut,tau_ratqs
1001!$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
1002      real zpt_conv(klon,klev)
1003
1004!c Parametres lies au nouveau schema de nuages (SB, PDF)
1005      real fact_cldcon
1006      real facttemps
1007      logical ok_newmicro
1008      save ok_newmicro
1009! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
1010!      real ref_liq(klon,klev), ref_ice(klon,klev)
1011
1012!$OMP THREADPRIVATE(ok_newmicro)
1013      save fact_cldcon,facttemps
1014!$OMP THREADPRIVATE(fact_cldcon,facttemps)
1015
1016      integer iflag_cldcon
1017      save iflag_cldcon
1018!$OMP THREADPRIVATE(iflag_cldcon)
1019! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
1020!      logical ptconv(klon,klev)
1021!IM cf. AM 081204 BEG
1022      logical ptconvth(klon,klev)
1023!IM cf. AM 081204 END
1024!c
1025!c Variables liees a l'ecriture de la bande histoire physique
1026!c
1027!c======================================================================
1028!c
1029!IM cf. AM 081204 BEG
1030!c   declarations pour sortir sur une sous-region
1031      integer imin_ins,imax_ins,jmin_ins,jmax_ins
1032      save imin_ins,imax_ins,jmin_ins,jmax_ins
1033!$OMP THREADPRIVATE(imin_ins,imax_ins,jmin_ins,jmax_ins)
1034!c      real lonmin_ins,lonmax_ins,latmin_ins
1035!c     s  ,latmax_ins
1036!c     data lonmin_ins,lonmax_ins,latmin_ins
1037!c    s  ,latmax_ins/
1038!c valeurs initiales     s   -5.,20.,41.,55./   
1039!c    s   100.,130.,-20.,20./
1040!c    s   -180.,180.,-90.,90./
1041!c======================================================================
1042!IM cf. AM 081204 END
1043
1044!c
1045      integer itau_w   ! pas de temps ecriture = itap + itau_phy
1046!c
1047!c
1048!c Variables locales pour effectuer les appels en serie
1049!c
1050! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule\n
1051!      REAL zx_rh(klon,klev)
1052!IM RH a 2m (la surface)
1053! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
1054!      REAL rh2m(klon), qsat2m(klon)
1055!      REAL tpot(klon), tpote(klon)
1056      REAL Lheat
1057
1058      INTEGER        length
1059      PARAMETER    ( length = 100 )
1060      REAL tabcntr0( length       )
1061!c
1062      INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev)
1063!IM
1064      INTEGER ndex2d1(iwmax)
1065!c
1066!IM AMIP2 BEG
1067      REAL moyglo, mountor
1068!IM 141004 BEG
1069      REAL zustrdr(klon), zvstrdr(klon)
1070      REAL zustrli(klon), zvstrli(klon)
1071      REAL zustrph(klon), zvstrph(klon)
1072      REAL zustrhi(klon), zvstrhi(klon)
1073      REAL aam, torsfc
1074!IM 141004 END
1075!IM 190504 BEG
1076      INTEGER ij, imp1jmp1
1077      PARAMETER(imp1jmp1=(iim+1)*jjmp1)
1078!cym A voir plus tard
1079      REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
1080      REAL padyn(iim+1,jjmp1,klev+1)
1081      REAL dudyn(iim+1,jjmp1,klev)
1082      REAL rlatdyn(iim+1,jjmp1)
1083!IM 190504 END
1084      LOGICAL ok_msk
1085      REAL msk(klon)
1086!IM
1087      REAL airetot, pi
1088!cym A voir plus tard
1089!cym      REAL zm_wo(jjmp1, klev)
1090!IM AMIP2 END
1091!c
1092      REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
1093      REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
1094      REAL zx_tmp_fi3d1(klon,klev+1) !variable temporaire pour champs 3D (kelvp1)
1095      REAL(KIND=8) zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
1096      REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
1097      REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
1098!c
1099      INTEGER nid_day, nid_mth, nid_ins, nid_mthnmc, nid_daynmc
1100      INTEGER nid_hfnmc, nid_day_seri, nid_ctesGCM
1101      SAVE nid_day, nid_mth, nid_ins, nid_mthnmc, nid_daynmc
1102      SAVE nid_hfnmc, nid_day_seri, nid_ctesGCM
1103!$OMP THREADPRIVATE(nid_day, nid_mth, nid_ins)
1104!$OMP THREADPRIVATE(nid_mthnmc, nid_daynmc, nid_hfnmc)
1105!$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
1106!c
1107!IM 280405 BEG
1108      INTEGER nid_bilKPins, nid_bilKPave
1109      SAVE nid_bilKPins, nid_bilKPave
1110!$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
1111!c
1112      REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
1113      REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
1114      REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
1115      REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
1116!c
1117      INTEGER nhori, nvert, nvert1, nvert3
1118      REAL zsto, zsto1, zsto2
1119      REAL zstophy, zstorad, zstohf, zstoday, zstomth, zout
1120      REAL zcals(napisccp), zcalh(napisccp), zoutj(napisccp)
1121      REAL zout_isccp(napisccp)
1122      SAVE zcals, zcalh, zoutj, zout_isccp
1123!$OMP THREADPRIVATE(zcals, zcalh, zoutj, zout_isccp)
1124
1125      real zjulian
1126      save zjulian
1127!$OMP THREADPRIVATE(zjulian)
1128
1129      character*20 modname
1130      character*80 abort_message
1131      logical ok_sync
1132      real date0
1133      integer idayref
1134
1135!C essai writephys
1136      integer fid_day, fid_mth, fid_ins
1137      parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
1138      integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
1139      parameter (prof2d_on = 1, prof3d_on = 2,                                       &
1140       &           prof2d_av = 3, prof3d_av = 4)
1141      character*30 nom_fichier
1142      character*40 varname
1143      character*40 vartitle
1144      character*20 varunits
1145!C     Variables liees au bilan d'energie et d'enthalpi
1146      REAL ztsol(klon)
1147      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot                           &
1148       &        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1149      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot                           &
1150       &        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1151!$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot,
1152!$OMP+              h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
1153      REAL      d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
1154      REAL      d_h_vcol_phy
1155      REAL      fs_bound, fq_bound
1156      SAVE      d_h_vcol_phy
1157!$OMP THREADPRIVATE(d_h_vcol_phy)
1158      REAL      zero_v(klon)
1159      CHARACTER*15 ztit
1160      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
1161      SAVE      ip_ebil
1162      DATA      ip_ebil/0/
1163!$OMP THREADPRIVATE(ip_ebil)
1164      INTEGER   if_ebil ! level for energy conserv. dignostics
1165      SAVE      if_ebil
1166!$OMP THREADPRIVATE(if_ebil)
1167!c+jld ec_conser
1168      REAL ZRCPD
1169!c-jld ec_conser
1170! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
1171!      REAL t2m(klon,nbsrf)  ! temperature a 2m
1172      REAL q2m(klon,nbsrf)  ! humidite a 2m
1173
1174!IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
1175! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
1176!      REAL zt2m(klon), zq2m(klon)             !temp., hum. 2m moyenne s/ 1 maille
1177!      REAL zustar(klon),zu10m(klon), zv10m(klon)  ! u* et vents a 10m moyennes s/1 maille
1178      CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
1179      CHARACTER*40 tinst, tave, typeval
1180      REAL cldtaupi(klon,klev)  ! Cloud optical thickness for pre-industrial (pi) aerosols
1181
1182! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
1183!      REAL re(klon, klev)       ! Cloud droplet effective radius
1184!      REAL fl(klon, klev)  ! denominator of re
1185
1186      REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
1187
1188      ! Aerosol optical properties
1189      CHARACTER*4, DIMENSION(naero_grp) :: rfname
1190      REAL, DIMENSION(klon)          :: aerindex     ! POLDER aerosol index
1191      REAL, DIMENSION(klon,klev)     :: mass_solu_aero    ! total mass concentration for all soluble aerosols[ug/m3]
1192      REAL, DIMENSION(klon,klev)     :: mass_solu_aero_pi ! - " - (pre-industrial value)
1193      INTEGER :: naero ! aerosol species
1194
1195      ! Parameters
1196      LOGICAL ok_ade, ok_aie    ! Apply aerosol (in)direct effects or not
1197      LOGICAL ok_cdnc          ! ok cloud droplet number concentration (O. Boucher 01-2013)
1198      REAL bl95_b0, bl95_b1   ! Parameter in Boucher and Lohmann (1995)
1199      SAVE ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1
1200!$OMP THREADPRIVATE(ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1)
1201      LOGICAL, SAVE :: aerosol_couple ! true  : calcul des aerosols dans INCA
1202                                      ! false : lecture des aerosol dans un fichier
1203!$OMP THREADPRIVATE(aerosol_couple)   
1204      INTEGER, SAVE :: flag_aerosol
1205!$OMP THREADPRIVATE(flag_aerosol)
1206      LOGICAL, SAVE :: new_aod
1207!$OMP THREADPRIVATE(new_aod)
1208!c
1209!c--STRAT AEROSOL
1210      LOGICAL, SAVE :: flag_aerosol_strat
1211!$OMP THREADPRIVATE(flag_aerosol_strat)
1212!cc-fin STRAT AEROSOL
1213!c
1214!c Declaration des constantes et des fonctions thermodynamiques
1215!c
1216      LOGICAL,SAVE :: first=.true.
1217!$OMP THREADPRIVATE(first)
1218
1219      integer iunit
1220
1221      integer, save::  read_climoz ! read ozone climatology
1222!C     (let it keep the default OpenMP shared attribute)
1223!C     Allowed values are 0, 1 and 2
1224!C     0: do not read an ozone climatology
1225!C     1: read a single ozone climatology that will be used day and night
1226!C     2: read two ozone climatologies, the average day and night
1227!C     climatology and the daylight climatology
1228
1229      integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
1230!C     (let it keep the default OpenMP shared attribute)
1231
1232      real, pointer, save:: press_climoz(:)
1233!C     (let it keep the default OpenMP shared attribute)
1234!     edges of pressure intervals for ozone climatologies, in Pa, in strictly
1235!     ascending order
1236
1237      integer, save:: co3i = 0
1238!     time index in NetCDF file of current ozone fields
1239!$OMP THREADPRIVATE(co3i)
1240
1241      integer ro3i
1242!     required time index in NetCDF file for the ozone fields, between 1
1243!     and 360
1244
1245      INTEGER ierr
1246#include "YOMCST.h"
1247#include "YOETHF.h"
1248#include "FCTTRE.h"
1249!IM 100106 BEG : pouvoir sortir les ctes de la physique
1250#include "conema3.h"
1251#include "fisrtilp.h"
1252#include "nuage.h"
1253#include "compbl.h"
1254!IM 100106 END : pouvoir sortir les ctes de la physique
1255!c
1256!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1257!c Declarations pour Simulateur COSP
1258!c============================================================
1259      real :: mr_ozone(klon,klev)
1260
1261!IM sorties fichier 1D paramLMDZ_phy.nc
1262      REAL :: zx_tmp_0d(1,1)
1263      INTEGER, PARAMETER :: np=1
1264      REAL,dimension(klon_glo)        :: rlat_glo
1265      REAL,dimension(klon_glo)        :: rlon_glo
1266      REAL gbils(1), gevap(1), gevapt(1), glat(1), gnet0(1), gnet(1)
1267      REAL grain(1), gtsol(1), gt2m(1), gprw(1)
1268
1269!IM stations CFMIP
1270      INTEGER, SAVE :: nCFMIP
1271!$OMP THREADPRIVATE(nCFMIP)
1272      INTEGER, PARAMETER :: npCFMIP=120
1273      INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
1274      REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
1275!$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP)
1276      INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
1277      REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
1278!$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
1279      INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
1280!$OMP THREADPRIVATE(iGCM, jGCM)
1281      logical, dimension(nfiles)            :: phys_out_filestations
1282      logical, parameter :: lNMC=.FALSE.
1283
1284!IM betaCRF
1285      REAL, SAVE :: pfree, beta_pbl, beta_free
1286!$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
1287      REAL, SAVE :: lon1_beta,  lon2_beta, lat1_beta, lat2_beta
1288!$OMP THREADPRIVATE(lon1_beta,  lon2_beta, lat1_beta, lat2_beta)
1289      LOGICAL, SAVE :: mskocean_beta
1290!$OMP THREADPRIVATE(mskocean_beta)
1291      REAL, dimension(klon, klev) :: beta         ! facteur sur cldtaurad et cldemirad pour evaluer les retros liees aux CRF
1292      REAL, dimension(klon, klev) :: cldtaurad    ! epaisseur optique pour radlwsw pour tester "CRF off"
1293      REAL, dimension(klon, klev) :: cldtaupirad  ! epaisseur optique pour radlwsw pour tester "CRF off"
1294      REAL, dimension(klon, klev) :: cldemirad    ! emissivite pour radlwsw pour tester "CRF off"
1295      REAL, dimension(klon, klev) :: cldfrarad    ! fraction nuageuse
1296
1297      INTEGER :: nbtr_tmp ! Number of tracer inside concvl
1298      REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
1299      integer iostat
1300
1301! L. Fita, LMD. Point for checkings...
1302      INTEGER                                            :: llp
1303      llp = 449
1304
1305!c======================================================================
1306! Gestion calendrier : mise a jour du module phys_cal_mod
1307!
1308      CALL phys_cal_update(jD_cur,jH_cur)
1309
1310!c======================================================================
1311! Ecriture eventuelle d'un profil verticale en entree de la physique.
1312! Utilise notamment en 1D mais peut etre active egalement en 3D
1313! en imposant la valeur de igout.
1314!c======================================================================
1315
1316      if (prt_level.ge.1) then
1317! Lluis
1318!          igout=klon/2+1/klon
1319          igout=llp
1320         write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1321         write(lunout,*)                                                             &
1322       & 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1323         write(lunout,*)                                                             &
1324       &  nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
1325
1326         write(lunout,*) 'paprs, play, phi, u, v, t'
1327         do k=1,klev
1328            write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),             &
1329       &   u(igout,k),v(igout,k),t(igout,k)
1330         enddo
1331         write(lunout,*) 'ovap (g/kg),  oliq (g/kg)'
1332         do k=1,klev
1333            write(lunout,*) qx(igout,k,1)*1000,qx(igout,k,2)*1000.
1334         enddo
1335      endif
1336
1337!c======================================================================
1338
1339!cym => necessaire pour iflag_con != 2   
1340      pmfd(:,:) = 0.
1341      pen_u(:,:) = 0.
1342      pen_d(:,:) = 0.
1343      pde_d(:,:) = 0.
1344      pde_u(:,:) = 0.
1345      aam=0.
1346
1347      torsfc=0.
1348
1349      forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
1350
1351      if (first) then
1352     
1353!cCR:nvelles variables convection/poches froides
1354     
1355      print*, '================================================='
1356      print*, 'Allocation des variables locales et sauvegardees'
1357
1358      call phys_local_var_init
1359
1360!c
1361      pasphys=pdtphys
1362!c     appel a la lecture du run.def physique
1363      call conf_phys(ok_journe, ok_mensuel,                                          &
1364       &     ok_instan, ok_hf,                                                       &
1365       &     ok_LES,                                                                 &
1366       &     callstats,                                                              &
1367       &     solarlong0,seuil_inversion,                                             &
1368       &     fact_cldcon, facttemps,ok_newmicro,iflag_radia,                         &
1369       &     iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,                  &
1370       &     ok_ade, ok_aie, ok_cdnc, aerosol_couple,                                &
1371       &     flag_aerosol, flag_aerosol_strat, new_aod,                              &
1372       &     bl95_b0, bl95_b1,                                                       &
1373!c     nv flags pour la convection et les poches froides                             &
1374       &     read_climoz,                                                            &
1375       &     alp_offset)
1376
1377
1378!!      call phys_state_var_init(read_climoz)
1379      call phys_output_var_init
1380
1381
1382!L. Fita, LMD. November 2013
1383! Initializing
1384! Variables have been already initialilzed from the WRF variables.
1385!   Using on purpose subroutines from physiq_limit_variables_mod
1386!!      CALL neige_initialize()
1387!!      CALL limit_initialize()
1388!!      CALL vars_limit_init(klon)
1389
1390! L. Fita, LMD January 2014. Initializing output variables
1391      CALL init_ovars_lmdz_NOmodule(klon,klev,nbsrf)
1392
1393! Allocating restart variables
1394!      ALLOCATE(qsol(klon))
1395!      ALLOCATE(fder(klon))
1396!      ALLOCATE(snow(klon, nbsrf))
1397!      ALLOCATE(qsurf(klon, nbsrf))
1398!      ALLOCATE(evap(klon, nbsrf))
1399!      ALLOCATE(rugos(klon, nbsrf))
1400!      ALLOCATE(agesno(klon, nbsrf))
1401!      ALLOCATE(ftsoil(klon, nsoilmx, nbsrf))
1402!!      ALLOCATE(restart_runoff(klon))
1403
1404!      CALL pbl_surface_init(qsol_rst, fder_rst, snow_rst, qsurf_rst, evap_rst,       &
1405!        rugos_rst, agesno_rst, ftsoil_rst)
1406!      CALL fonte_neige_init(restart_runoff)
1407!      qsol=qsol_rst
1408
1409      print*, '================================================='
1410!c
1411          dnwd0=0.0
1412          ftd=0.0
1413          fqd=0.0
1414          cin=0.
1415!cym Attention pbase pas initialise dans concvl !!!!
1416          pbase=0
1417!IM 180608
1418
1419        itau_con=0
1420        first=.false.
1421
1422      endif  ! first
1423
1424       modname = 'physiq'
1425!IM
1426      IF (ip_ebil_phy.ge.1) THEN
1427        DO i=1,klon
1428          zero_v(i)=0.
1429        END DO
1430      END IF
1431      ok_sync=.TRUE.
1432
1433      IF (debut) THEN
1434         CALL suphel ! initialiser constantes et parametres phys.
1435      ENDIF
1436
1437      if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
1438
1439
1440!c======================================================================
1441! Gestion calendrier : mise a jour du module phys_cal_mod
1442!
1443!c     CALL phys_cal_update(jD_cur,jH_cur)
1444
1445!c
1446!c Si c'est le debut, il faut initialiser plusieurs choses
1447!c          ********
1448!c
1449
1450       IF (debut) THEN
1451!rv
1452!cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
1453!cde la convection a partir des caracteristiques du thermique
1454         wght_th(:,:)=1.
1455         lalim_conv(:)=1
1456!cRC
1457         ustar(:,:)=0.
1458         u10m(:,:)=0.
1459         v10m(:,:)=0.
1460         rain_con(:)=0.
1461         snow_con(:)=0.
1462         topswai(:)=0.
1463         topswad(:)=0.
1464         solswai(:)=0.
1465         solswad(:)=0.
1466
1467         wmax_th(:)=0.
1468         tau_overturning_th(:)=0.
1469
1470         IF (type_trac == 'inca') THEN
1471            ! jg : initialisation jusqu'au ces variables sont dans restart
1472            ccm(:,:,:) = 0.
1473            tau_aero(:,:,:,:) = 0.
1474            piz_aero(:,:,:,:) = 0.
1475            cg_aero(:,:,:,:) = 0.
1476         END IF
1477
1478         rnebcon0(:,:) = 0.0
1479         clwcon0(:,:) = 0.0
1480         rnebcon(:,:) = 0.0
1481         clwcon(:,:) = 0.0
1482
1483!IM     
1484         IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1485!c
1486      print*,'iflag_coupl,iflag_clos,iflag_wake',                                    &
1487       &   iflag_coupl,iflag_clos,iflag_wake
1488      print*,'CYCLE_DIURNE', cycle_diurne
1489!c
1490      IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
1491         abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
1492         CALL abort_gcm (modname,abort_message,1)
1493      ENDIF
1494!c
1495      IF(ok_isccp.AND.iflag_con.LE.2) THEN
1496         abort_message = 'ISCCP-like outputs may be available for KE                 &
1497       &(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
1498         CALL abort_gcm (modname,abort_message,1)
1499      ENDIF
1500!c
1501!c Initialiser les compteurs:
1502!c
1503         itap    = 0
1504         itaprad = 0
1505
1506!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1507!! Un petit travail \`a faire ici.
1508!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1509
1510         if (iflag_pbl>1) then
1511             PRINT*, "Using method MELLOR&YAMADA"
1512         endif
1513
1514!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1515! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1516! dyn3d
1517! Attention : la version precedente n'etait pas tres propre.
1518! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1519! pour obtenir le meme resultat.
1520         dtime=pdtphys
1521         radpas = NINT( 86400./dtime/nbapp_rad)
1522!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1523
1524         CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1525
1526         IF (klon_glo==1) THEN
1527         coefh=0. ; coefm=0. ; pbl_tke=0.
1528         coefh(:,2,:)=1.e-2 ; coefm(:,2,:)=1.e-2 ; pbl_tke(:,2,:)=1.e-2
1529         PRINT*,'FH WARNING : lignes a supprimer'
1530         ENDIF
1531!IM begin
1532          print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)             &
1533       &,ratqs(1,1)
1534!IM end
1535
1536
1537
1538!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1539!c
1540!C on remet le calendrier a zero
1541!c
1542         IF (raz_date .eq. 1) THEN
1543           itau_phy = 0
1544         ENDIF
1545
1546!IM cf. AM 081204 BEG
1547         PRINT*,'cycle_diurne3 =',cycle_diurne
1548!IM cf. AM 081204 END
1549!c
1550         CALL printflag( tabcntr0,radpas,ok_journe,                                  &
1551       &                    ok_instan, ok_region )
1552!c
1553         IF (ABS(dtime-pdtphys).GT.0.001) THEN
1554            WRITE(lunout,*) 'Pas physique n est pas correct',dtime,                  &
1555       &                        pdtphys
1556            abort_message='Pas physique n est pas correct '
1557!           call abort_gcm(modname,abort_message,1)
1558            dtime=pdtphys
1559         ENDIF
1560         IF (nlon .NE. klon) THEN
1561            WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,               &
1562       &                      klon
1563            abort_message='nlon et klon ne sont pas coherents'
1564            call abort_gcm(modname,abort_message,1)
1565         ENDIF
1566         IF (nlev .NE. klev) THEN
1567            WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,               &
1568       &                       klev
1569            abort_message='nlev et klev ne sont pas coherents'
1570            call abort_gcm(modname,abort_message,1)
1571         ENDIF
1572!c
1573         IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
1574           WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1575           WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1576           abort_message='Nbre d appels au rayonnement insuffisant'
1577           call abort_gcm(modname,abort_message,1)
1578         ENDIF
1579         WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1580         WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",             &
1581       &                   ok_cvl
1582!c
1583!cKE43
1584!c Initialisation pour la convection de K.E. (sb):
1585         IF (iflag_con.GE.3) THEN
1586
1587         WRITE(lunout,*)"*** Convection de Kerry Emanuel 4.3  "
1588         WRITE(lunout,*)                                                             &
1589       &      "On va utiliser le melange convectif des traceurs qui"
1590         WRITE(lunout,*)"est calcule dans convect4.3"
1591         WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1592
1593          DO i = 1, klon
1594           ema_cbmf(i) = 0.
1595           ema_pcb(i)  = 0.
1596           ema_pct(i)  = 0.
1597!c          ema_workcbmf(i) = 0.
1598          ENDDO
1599!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1600          DO i = 1, klon
1601           ibas_con(i) = 1
1602           itop_con(i) = 1
1603          ENDDO
1604!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1605!c===============================================================================
1606!cCR:04.12.07: initialisations poches froides
1607!c Controle de ALE et ALP pour la fermeture convective (jyg)
1608          if (iflag_wake>=1) then
1609            CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr                &
1610       &                  ,alp_bl_prescr, ale_bl_prescr)
1611!c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1612!c        print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1613          endif
1614
1615        do i = 1,klon
1616         Ale_bl(i)=0.
1617         Alp_bl(i)=0.
1618        enddo
1619
1620!c================================================================================
1621!IM stations CFMIP
1622      nCFMIP=npCFMIP
1623      OPEN(98,file='npCFMIP_param.data',status='old',                                &
1624       &          form='formatted',iostat=iostat)
1625            if (iostat == 0) then
1626      READ(98,*,end=998) nCFMIP
1627998   CONTINUE
1628      CLOSE(98)
1629      CONTINUE
1630      IF(nCFMIP.GT.npCFMIP) THEN
1631       print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
1632       CALL abort
1633      else
1634       print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
1635      ENDIF
1636
1637!c
1638      ALLOCATE(tabCFMIP(nCFMIP))
1639      ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
1640      ALLOCATE(tabijGCM(nCFMIP))
1641      ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
1642      ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
1643!c
1644!c lecture des nCFMIP stations CFMIP, de leur numero
1645!c et des coordonnees geographiques lonCFMIP, latCFMIP
1646!c
1647         CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP,                           &
1648       &lonCFMIP, latCFMIP)
1649!c
1650!c identification des
1651!c 1) coordonnees lonGCM, latGCM des points CFMIP dans la grille de LMDZ
1652!c 2) indices points tabijGCM de la grille physique 1d sur klon points
1653!c 3) indices iGCM, jGCM de la grille physique 2d
1654!c
1655         CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP,                 &
1656       &tabijGCM, lonGCM, latGCM, iGCM, jGCM)
1657!c
1658            else
1659               ALLOCATE(tabijGCM(0))
1660               ALLOCATE(lonGCM(0), latGCM(0))
1661               ALLOCATE(iGCM(0), jGCM(0))
1662            end if
1663         else
1664            ALLOCATE(tabijGCM(0))
1665            ALLOCATE(lonGCM(0), latGCM(0))
1666            ALLOCATE(iGCM(0), jGCM(0))
1667         ENDIF
1668 
1669           DO i=1,klon
1670             rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1671           ENDDO
1672
1673!c34EK
1674         IF (ok_orodr) THEN
1675
1676!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1677! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1678! justement quand ok_orodr = false.
1679! ce rugoro est utilise par la couche limite et fait double emploi
1680! avec les param\'etrisations sp\'ecifiques de Francois Lott.
1681!           DO i=1,klon
1682!             rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1683!           ENDDO
1684!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1685           IF (ok_strato) THEN
1686             CALL SUGWD_strato(klon,klev,paprs,pplay)
1687           ELSE
1688             CALL SUGWD(klon,klev,paprs,pplay)
1689           ENDIF
1690           
1691           DO i=1,klon
1692             zuthe(i)=0.
1693             zvthe(i)=0.
1694             if(zstd(i).gt.10.)then
1695               zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1696               zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1697             endif
1698           ENDDO
1699         ENDIF
1700!c
1701!c
1702         lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
1703         WRITE(lunout,*)'La frequence de lecture surface est de ',                   &
1704       &                   lmt_pas
1705!c
1706      capemaxcels = 't_max(X)'
1707      t2mincels = 't_min(X)'
1708      t2maxcels = 't_max(X)'
1709      tinst = 'inst(X)'
1710      tave = 'ave(X)'
1711!IM cf. AM 081204 BEG
1712      write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1713!IM cf. AM 081204 END
1714!c
1715!c=============================================================
1716!c   Initialisation des sorties
1717!c=============================================================
1718
1719#ifdef CPP_IOIPSL
1720
1721!$OMP MASTER
1722       call phys_output_open(rlon,rlat,nCFMIP,tabijGCM,                              &
1723       &                       iGCM,jGCM,lonGCM,latGCM,                              &
1724       &                       jjmp1,nlevSTD,clevSTD,                                &
1725       &                       nbteta, ctetaSTD, dtime,ok_veget,                     &
1726       &                       type_ocean,iflag_pbl,ok_mensuel,ok_journe,            &
1727       &                       ok_hf,ok_instan,ok_LES,ok_ade,ok_aie,                 &
1728       &                       read_climoz, phys_out_filestations,                   &
1729       &                       new_aod, aerosol_couple,                              &
1730       &                       flag_aerosol_strat )
1731!$OMP END MASTER
1732!$OMP BARRIER
1733
1734#ifdef histISCCP
1735#include "ini_histISCCP.h"
1736#endif
1737
1738#ifdef histNMC
1739#include "ini_histhfNMC.h"
1740#include "ini_histdayNMC.h"
1741#include "ini_histmthNMC.h"
1742#endif
1743
1744#include "ini_histday_seri.h"
1745
1746#include "ini_paramLMDZ_phy.h"
1747
1748#endif
1749         ecrit_reg = ecrit_reg * un_jour
1750         ecrit_tra = ecrit_tra * un_jour
1751       
1752!cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1753      date0 = jD_ref
1754      WRITE(*,*) 'physiq date0 : ',date0
1755!c
1756!c
1757!c
1758!c Prescrire l'ozone dans l'atmosphere
1759!c
1760!c
1761!cc         DO i = 1, klon
1762!cc         DO k = 1, klev
1763!cc            CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1764!cc         ENDDO
1765!cc         ENDDO
1766!c
1767      IF (type_trac == 'inca') THEN
1768#ifdef INCA
1769         CALL VTe(VTphysiq)
1770         CALL VTb(VTinca)
1771!         iii = MOD(NINT(xjour),360)
1772!         calday = REAL(iii) + jH_cur
1773         calday = REAL(days_elapsed) + jH_cur
1774         WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1775
1776         CALL chemini(                                                               &
1777       &                   rg,                                                       &
1778       &                   ra,                                                       &
1779       &                   airephy,                                                  &
1780       &                   rlat,                                                     &
1781       &                   rlon,                                                     &
1782       &                   presnivs,                                                 &
1783       &                   calday,                                                   &
1784       &                   klon,                                                     &
1785       &                   nqtot,                                                    &
1786       &                   pdtphys,                                                  &
1787       &                   annee_ref,                                                &
1788       &                   day_ref,                                                  &
1789       &                   itau_phy)
1790
1791         CALL VTe(VTinca)
1792         CALL VTb(VTphysiq)
1793#endif
1794      END IF
1795!c
1796!c
1797!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1798! Nouvelle initialisation pour le rayonnement RRTM
1799!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1800
1801      call iniradia(klon,klev,paprs(1,1:klev+1))
1802
1803!$OMP single
1804      if (read_climoz >= 1) then
1805         call open_climoz(ncid_climoz, press_climoz)
1806      END IF
1807!$OMP end single
1808!c
1809!IM betaCRF
1810      pfree=70000. !Pa
1811      beta_pbl=1.
1812      beta_free=1.
1813      lon1_beta=-180.
1814      lon2_beta=+180.
1815      lat1_beta=90.
1816      lat2_beta=-90.
1817      mskocean_beta=.FALSE.
1818
1819      OPEN(99,file='beta_crf.data',status='old',                                     &
1820       &          form='formatted',err=9999)
1821      READ(99,*,end=9998) pfree
1822      READ(99,*,end=9998) beta_pbl
1823      READ(99,*,end=9998) beta_free
1824      READ(99,*,end=9998) lon1_beta
1825      READ(99,*,end=9998) lon2_beta
1826      READ(99,*,end=9998) lat1_beta
1827      READ(99,*,end=9998) lat2_beta
1828      READ(99,*,end=9998) mskocean_beta
18299998  Continue
1830      CLOSE(99)
18319999  Continue
1832      WRITE(*,*)'pfree=',pfree
1833      WRITE(*,*)'beta_pbl=',beta_pbl
1834      WRITE(*,*)'beta_free=',beta_free
1835      WRITE(*,*)'lon1_beta=',lon1_beta
1836      WRITE(*,*)'lon2_beta=',lon2_beta
1837      WRITE(*,*)'lat1_beta=',lat1_beta
1838      WRITE(*,*)'lat2_beta=',lat2_beta
1839      WRITE(*,*)'mskocean_beta=',mskocean_beta
1840      ENDIF
1841!
1842!   ****************     Fin  de   IF ( debut  )   ***************
1843!
1844!
1845! Incrementer le compteur de la physique
1846!
1847      itap   = itap + 1
1848!c
1849!
1850! Update fraction of the sub-surfaces (pctsrf) and
1851! initialize, where a new fraction has appeared, all variables depending
1852! on the surface fraction.
1853!
1854      CALL change_srf_frac(itap, dtime, days_elapsed+1,                              &
1855       &     pctsrf, falb1, falb2, ftsol, ustar, u10m, v10m, pbl_tke)
1856
1857! Update time and other variables in Reprobus
1858      IF (type_trac == 'repr') THEN
1859#ifdef REPROBUS
1860         CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
1861         print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
1862         CALL Rtime(debut)
1863#endif
1864      END IF
1865
1866
1867! Tendances bidons pour les processus qui n'affectent pas certaines
1868! variables.
1869      du0(:,:)=0.
1870      dv0(:,:)=0.
1871      dq0(:,:)=0.
1872      dql0(:,:)=0.
1873!c
1874!c Mettre a zero des variables de sortie (pour securite)
1875!c
1876      DO i = 1, klon
1877         d_ps(i) = 0.0
1878      ENDDO
1879      DO k = 1, klev
1880      DO i = 1, klon
1881         d_t(i,k) = 0.0
1882         d_u(i,k) = 0.0
1883         d_v(i,k) = 0.0
1884      ENDDO
1885      ENDDO
1886      DO iq = 1, nqtot
1887      DO k = 1, klev
1888      DO i = 1, klon
1889         d_qx(i,k,iq) = 0.0
1890      ENDDO
1891      ENDDO
1892      ENDDO
1893      da(:,:)=0.
1894      mp(:,:)=0.
1895      phi(:,:,:)=0.
1896! RomP >>>
1897      phi2(:,:,:)=0.
1898      beta_prec_fisrt(:,:)=0.
1899      beta_prec(:,:)=0.
1900      epmlmMm(:,:,:)=0.
1901      eplaMm(:,:)=0.
1902      d1a(:,:)=0.
1903      dam(:,:)=0.
1904          pmflxr=0.
1905          pmflxs=0.
1906! RomP <<<
1907
1908!c
1909!c Ne pas affecter les valeurs entrees de u, v, h, et q
1910!c
1911      DO k = 1, klev
1912      DO i = 1, klon
1913         t_seri(i,k)  = t(i,k)
1914         u_seri(i,k)  = u(i,k)
1915         v_seri(i,k)  = v(i,k)
1916         q_seri(i,k)  = qx(i,k,ivap)
1917         ql_seri(i,k) = qx(i,k,iliq)
1918         qs_seri(i,k) = 0.
1919      ENDDO
1920      ENDDO
1921      tke0(:,:)=pbl_tke(:,:,is_ave)
1922      IF (nqtot.GE.3) THEN
1923      DO iq = 3, nqtot
1924      DO  k = 1, klev
1925      DO  i = 1, klon
1926         tr_seri(i,k,iq-2) = qx(i,k,iq)
1927      ENDDO
1928      ENDDO
1929      ENDDO
1930      ELSE
1931      DO k = 1, klev
1932      DO i = 1, klon
1933         tr_seri(i,k,1) = 0.0
1934      ENDDO
1935      ENDDO
1936      ENDIF
1937!C
1938      DO i = 1, klon
1939        ztsol(i) = 0.
1940      ENDDO
1941      DO nsrf = 1, nbsrf
1942        DO i = 1, klon
1943          ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1944        ENDDO
1945      ENDDO
1946!IM
1947      IF (ip_ebil_phy.ge.1) THEN
1948        ztit='after dynamic'
1949        CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime                             &
1950       &      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay              &
1951       &      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1952!C     Comme les tendances de la physique sont ajoute dans la dynamique,
1953!C     on devrait avoir que la variation d'entalpie par la dynamique
1954!C     est egale a la variation de la physique au pas de temps precedent.
1955!C     Donc la somme de ces 2 variations devrait etre nulle.
1956        call diagphy(airephy,ztit,ip_ebil_phy                                        &
1957       &      , zero_v, zero_v, zero_v, zero_v, zero_v                               &
1958       &      , zero_v, zero_v, zero_v, ztsol                                        &
1959       &      , d_h_vcol+d_h_vcol_phy, d_qt, 0.                                      &
1960       &      , fs_bound, fq_bound )
1961      END IF
1962      PRINT *,'  Lluis 1 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
1963
1964!c Diagnostiquer la tendance dynamique
1965!c
1966      IF (ancien_ok) THEN
1967         DO k = 1, klev
1968         DO i = 1, klon
1969            d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1970            d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1971            d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1972            d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1973         ENDDO
1974         ENDDO
1975!!! RomP >>>   td dyn traceur
1976       IF (nqtot.GE.3) THEN
1977          DO iq = 3, nqtot
1978          DO k = 1, klev
1979          DO i = 1, klon
1980            d_tr_dyn(i,k,iq-2)=                                                      &
1981       &       (tr_seri(i,k,iq-2)-tr_ancien(i,k,iq-2))/dtime
1982!         iiq=niadv(iq)
1983!         print*,i,k," d_tr_dyn",d_tr_dyn(i,k,iq-2),"tra:",iq,tname(iiq)
1984          ENDDO
1985          ENDDO
1986          ENDDO
1987       ENDIF
1988!!! RomP <<<
1989      ELSE
1990         DO k = 1, klev
1991         DO i = 1, klon
1992            d_u_dyn(i,k) = 0.0
1993            d_v_dyn(i,k) = 0.0
1994            d_t_dyn(i,k) = 0.0
1995            d_q_dyn(i,k) = 0.0
1996         ENDDO
1997         ENDDO
1998!!! RomP >>>   td dyn traceur
1999        IF (nqtot.GE.3) THEN
2000          DO iq = 3, nqtot
2001          DO k = 1, klev
2002          DO i = 1, klon
2003            d_tr_dyn(i,k,iq-2)= 0.0
2004          ENDDO
2005          ENDDO
2006          ENDDO
2007       ENDIF
2008!!! RomP <<<
2009         ancien_ok = .TRUE.
2010      ENDIF
2011!c
2012!c Ajouter le geopotentiel du sol:
2013!c
2014      DO k = 1, klev
2015      DO i = 1, klon
2016         zphi(i,k) = pphi(i,k) + pphis(i)
2017      ENDDO
2018      ENDDO
2019!c
2020!c Verifier les temperatures
2021!c
2022!IM BEG
2023      IF (check) THEN
2024       amn=MIN(ftsol(1,is_ter),1000.)
2025       amx=MAX(ftsol(1,is_ter),-1000.)
2026       DO i=2, klon
2027        amn=MIN(ftsol(i,is_ter),amn)
2028        amx=MAX(ftsol(i,is_ter),amx)
2029       ENDDO
2030!c
2031       PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
2032      ENDIF !(check) THEN
2033
2034!IM END
2035!c
2036
2037      CALL hgardfou(t_seri,ftsol,'debutphy')
2038!c
2039!IM BEG
2040      IF (check) THEN
2041       amn=MIN(ftsol(1,is_ter),1000.)
2042       amx=MAX(ftsol(1,is_ter),-1000.)
2043       DO i=2, klon
2044        amn=MIN(ftsol(i,is_ter),amn)
2045        amx=MAX(ftsol(i,is_ter),amx)
2046       ENDDO
2047!c
2048       PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
2049      ENDIF !(check) THEN
2050
2051!IM END
2052!c
2053!c Mettre en action les conditions aux limites (albedo, sst, etc.).
2054!c Prescrire l'ozone et calculer l'albedo sur l'ocean.
2055!c
2056      if (read_climoz >= 1) then
2057!C        Ozone from a file
2058!        Update required ozone index:
2059         ro3i = int((days_elapsed + jh_cur - jh_1jan)                                &
2060       &        / ioget_year_len(year_cur) * 360.) + 1
2061         if (ro3i == 361) ro3i = 360
2062!C        (This should never occur, except perhaps because of roundup
2063!C        error. See documentation.)
2064         if (ro3i /= co3i) then
2065!C           Update ozone field:
2066            if (read_climoz == 1) then
2067               call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i,                 &
2068       &              press_in_edg=press_climoz, paprs=paprs, v3=wo)
2069            else
2070!C              read_climoz == 2
2071               call regr_pr_av(ncid_climoz,                                          &
2072       &              (/"tro3         ", "tro3_daylight"/),                          &
2073       &              julien=ro3i, press_in_edg=press_climoz, paprs=paprs,           &
2074       &              v3=wo)
2075            end if
2076!           Convert from mole fraction of ozone to column density of ozone in a
2077!           cell, in kDU:
2078            forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l)                    &
2079       &           * rmo3 / rmd * zmasse / dobson_u / 1e3
2080!C           (By regridding ozone values for LMDZ only once every 360th of
2081!C           year, we have already neglected the variation of pressure in one
2082!C           360th of year. So do not recompute "wo" at each time step even if
2083!C           "zmasse" changes a little.)
2084            co3i = ro3i
2085         end if
2086      elseif (MOD(itap-1,lmt_pas) == 0) THEN
2087!C        Once per day, update ozone from Royer:
2088         wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
2089      ENDIF
2090
2091!c
2092!c Re-evaporer l'eau liquide nuageuse
2093!c
2094      DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
2095      DO i = 1, klon
2096         zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
2097!c        zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
2098         zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
2099         zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
2100         zb = MAX(0.0,ql_seri(i,k))
2101         za = - MAX(0.0,ql_seri(i,k))                                                &
2102       &                  * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
2103         t_seri(i,k) = t_seri(i,k) + za
2104         q_seri(i,k) = q_seri(i,k) + zb
2105         ql_seri(i,k) = 0.0
2106         d_t_eva(i,k) = za
2107         d_q_eva(i,k) = zb
2108      ENDDO
2109      ENDDO
2110!IM
2111      IF (ip_ebil_phy.ge.2) THEN
2112        ztit='after reevap'
2113        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime                             &
2114       &      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay              &
2115       &      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2116         call diagphy(airephy,ztit,ip_ebil_phy                                       &
2117       &      , zero_v, zero_v, zero_v, zero_v, zero_v                               &
2118       &      , zero_v, zero_v, zero_v, ztsol                                        &
2119       &      , d_h_vcol, d_qt, d_ec                                                 &
2120       &      , fs_bound, fq_bound )
2121!C
2122      END IF
2123      PRINT *,'  Lluis 2 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
2124!c
2125!c=========================================================================
2126! Calculs de l'orbite.
2127! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2128! doit donc etre plac\'e avant radlwsw et pbl_surface
2129
2130!!!   jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2131      call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
2132      day_since_equinox = (jD_cur + jH_cur) - jD_eq
2133!
2134!   choix entre calcul de la longitude solaire vraie ou valeur fixee a
2135!   solarlong0
2136      if (solarlong0<-999.) then
2137       if (new_orbit) then
2138! calcul selon la routine utilisee pour les planetes
2139        call solarlong(day_since_equinox, zlongi, dist)
2140       else
2141! calcul selon la routine utilisee pour l'AR4
2142        CALL orbite(REAL(days_elapsed+1),zlongi,dist)
2143       endif
2144      else
2145           zlongi=solarlong0  ! longitude solaire vraie
2146           dist=1.            ! distance au soleil / moyenne
2147      endif
2148      if(prt_level.ge.1)                                                &
2149       &    write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist,&
2150       ' jD_cur: ',jD_cur,' jH_cur: ',jH_cur,' days_elapsed: ', &
2151       REAL(days_elapsed+1)
2152
2153
2154!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2155! Calcul de l'ensoleillement :
2156! ============================
2157! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2158! l'annee a partir d'une formule analytique.
2159! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
2160! non nul aux poles.
2161      IF (abs(solarlong0-1000.)<1.e-4) then
2162         call zenang_an(cycle_diurne,jH_cur,rlat,rlon,rmu0,fract)
2163      ELSE
2164!  Avec ou sans cycle diurne
2165         IF (cycle_diurne) THEN
2166           zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
2167           CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
2168         ELSE
2169           CALL angle(zlongi, rlat, fract, rmu0)
2170         ENDIF
2171      ENDIF
2172
2173      if (mydebug) then
2174        call writefield_phy('u_seri',u_seri,llm)
2175        call writefield_phy('v_seri',v_seri,llm)
2176        call writefield_phy('t_seri',t_seri,llm)
2177        call writefield_phy('q_seri',q_seri,llm)
2178      endif
2179
2180!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2181!c Appel au pbl_surface : Planetary Boudary Layer et Surface
2182!c Cela implique tous les interactions des sous-surfaces et la partie diffusion
2183!c turbulent du couche limit.
2184!c
2185!c Certains varibales de sorties de pbl_surface sont utiliser que pour
2186!c ecriture des fihiers hist_XXXX.nc, ces sont :
2187!c   qsol,      zq2m,      s_pblh,  s_lcl,
2188!c   s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
2189!c   s_therm,   s_trmb1,   s_trmb2, s_trmb3,
2190!c   zxrugs,    zu10m,     zv10m,   fder,
2191!c   zxqsurf,   rh2m,      zxfluxu, zxfluxv,
2192!c   frugs,     agesno,    fsollw,  fsolsw,
2193!c   d_ts,      fevap,     fluxlat, t2m,
2194!c   wfbils,    wfbilo,    fluxt,   fluxu, fluxv,
2195!c
2196!c Certains ne sont pas utiliser du tout :
2197!c   dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2198!c
2199
2200!c Calcul de l'humidite de saturation au niveau du sol
2201
2202
2203
2204      if (iflag_pbl/=0) then
2205
2206      CALL pbl_surface(                                                              &
2207       &     dtime,     date0,     itap,    days_elapsed+1,                          &
2208       &     debut,     lafin,                                                       &
2209       &     rlon,      rlat,      rugoro,  rmu0,                                    &
2210       &     rain_fall, snow_fall, solsw,   sollw,                                   &
2211       &     t_seri,    q_seri,    u_seri,  v_seri,                                  &
2212       &     pplay,     paprs,     pctsrf,                                           &
2213       &     ftsol,     falb1,     falb2,   ustar, u10m,   v10m,                     &
2214       &     sollwdown, cdragh,    cdragm,  u1,    v1,                               &
2215       &     albsol1,   albsol2,   sens,    evap,                                    &
2216       &     zxtsol,    zxfluxlat, zt2m,    qsat2m,                                  &
2217       &     d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf, d_t_diss,                       &
2218       &     coefh,     coefm,     slab_wfbils,                                      &
2219       &     qsol,      zq2m,      s_pblh,  s_lcl,                                   &
2220       &     s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,                                  &
2221       &     s_therm,   s_trmb1,   s_trmb2, s_trmb3,                                 &
2222       &     zxrugs,    zustar, zu10m,     zv10m,   fder,                            &
2223       &     zxqsurf,   rh2m,      zxfluxu, zxfluxv,                                 &
2224       &     frugs,     agesno,    fsollw,  fsolsw,                                  &
2225       &     d_ts,      fevap,     fluxlat, t2m,                                     &
2226       &     wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,                           &
2227       &     dsens,     devap,     zxsnow,                                           &
2228       &     zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
2229
2230!-----------------------------------------------------------------------------------------
2231! ajout des tendances de la diffusion turbulente
2232      CALL add_phys_tend                                                             &
2233       &     (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,'vdf')
2234!-----------------------------------------------------------------------------------------
2235      if (mydebug) then
2236        call writefield_phy('u_seri',u_seri,llm)
2237        call writefield_phy('v_seri',v_seri,llm)
2238        call writefield_phy('t_seri',t_seri,llm)
2239        call writefield_phy('q_seri',q_seri,llm)
2240      endif
2241
2242         CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh,                            &
2243       &      t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
2244
2245
2246      IF (ip_ebil_phy.ge.2) THEN
2247        ztit='after surface_main'
2248        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime                             &
2249       &      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay              &
2250       &      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2251         call diagphy(airephy,ztit,ip_ebil_phy                                       &
2252       &      , zero_v, zero_v, zero_v, zero_v, sens                                 &
2253       &      , evap  , zero_v, zero_v, ztsol                                        &
2254       &      , d_h_vcol, d_qt, d_ec                                                 &
2255       &      , fs_bound, fq_bound )
2256      END IF
2257
2258      ENDIF
2259      PRINT *,'  Lluis 3 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
2260!c =================================================================== c
2261!c   Calcul de Qsat
2262
2263      DO k = 1, klev
2264      DO i = 1, klon
2265         zx_t = t_seri(i,k)
2266         IF (thermcep) THEN
2267            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2268            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2269            zx_qs  = MIN(0.5,zx_qs)
2270            zcor   = 1./(1.-retv*zx_qs)
2271            zx_qs  = zx_qs*zcor
2272         ELSE
2273           IF (zx_t.LT.t_coup) THEN
2274              zx_qs = qsats(zx_t)/pplay(i,k)
2275           ELSE
2276              zx_qs = qsatl(zx_t)/pplay(i,k)
2277           ENDIF
2278         ENDIF
2279         zqsat(i,k)=zx_qs
2280      ENDDO
2281      ENDDO
2282
2283      if (prt_level.ge.1) then
2284      write(lunout,*) 'L   qsat (g/kg) avant clouds_gno'
2285      write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
2286      endif
2287!c
2288!c Appeler la convection (au choix)
2289!c
2290      DO k = 1, klev
2291      DO i = 1, klon
2292         conv_q(i,k) = d_q_dyn(i,k)                                                  &
2293       &               + d_q_vdf(i,k)/dtime
2294         conv_t(i,k) = d_t_dyn(i,k)                                                  &
2295       &               + d_t_vdf(i,k)/dtime
2296      ENDDO
2297      ENDDO
2298      IF (check) THEN
2299         za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2300         WRITE(lunout,*) "avantcon=", za
2301      ENDIF
2302      zx_ajustq = .FALSE.
2303      IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2304      IF (zx_ajustq) THEN
2305         DO i = 1, klon
2306            z_avant(i) = 0.0
2307         ENDDO
2308         DO k = 1, klev
2309         DO i = 1, klon
2310            z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))                     &
2311       &                        *(paprs(i,k)-paprs(i,k+1))/RG
2312         ENDDO
2313         ENDDO
2314      ENDIF
2315
2316!c Calcule de vitesse verticale a partir de flux de masse verticale
2317      DO k = 1, klev
2318         DO i = 1, klon
2319            omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
2320         END DO
2321      END DO
2322      if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ',                      &
2323       &     omega(igout, :)
2324
2325      IF (iflag_con.EQ.1) THEN
2326        abort_message ='reactiver le call conlmd dans physiq.F'
2327        CALL abort_gcm (modname,abort_message,1)
2328!c     CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2329!c    .             d_t_con, d_q_con,
2330!c    .             rain_con, snow_con, ibas_con, itop_con)
2331      ELSE IF (iflag_con.EQ.2) THEN
2332      CALL conflx(dtime, paprs, pplay, t_seri, q_seri,                               &
2333       &            conv_t, conv_q, -evap, omega,                                    &
2334       &            d_t_con, d_q_con, rain_con, snow_con,                            &
2335       &            pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,                          &
2336       &            kcbot, kctop, kdtop, pmflxr, pmflxs)
2337      d_u_con = 0.
2338      d_v_con = 0.
2339
2340      WHERE (rain_con < 0.) rain_con = 0.
2341      WHERE (snow_con < 0.) snow_con = 0.
2342      DO i = 1, klon
2343         ibas_con(i) = klev+1 - kcbot(i)
2344         itop_con(i) = klev+1 - kctop(i)
2345      ENDDO
2346      ELSE IF (iflag_con.GE.3) THEN
2347!c nb of tracers for the KE convection:
2348!c MAF la partie traceurs est faite dans phytrac
2349!c on met ntra=1 pour limiter les appels mais on peut
2350!c supprimer les calculs / ftra.
2351              ntra = 1
2352
2353!c=====================================================================================
2354!cajout pour la parametrisation des poches froides:
2355!ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
2356      do k=1,klev
2357            do i=1,klon
2358             if (iflag_wake>=1) then
2359             t_wake(i,k) = t_seri(i,k)                                               &
2360       &           +(1-wake_s(i))*wake_deltat(i,k)
2361             q_wake(i,k) = q_seri(i,k)                                               &
2362       &           +(1-wake_s(i))*wake_deltaq(i,k)
2363             t_undi(i,k) = t_seri(i,k)                                               &
2364       &           -wake_s(i)*wake_deltat(i,k)
2365             q_undi(i,k) = q_seri(i,k)                                               &
2366       &           -wake_s(i)*wake_deltaq(i,k)
2367             else
2368             t_wake(i,k) = t_seri(i,k)
2369             q_wake(i,k) = q_seri(i,k)
2370             t_undi(i,k) = t_seri(i,k)
2371             q_undi(i,k) = q_seri(i,k)
2372             endif
2373            enddo
2374         enddo
2375     
2376!cc--   Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
2377!cc--    pour le soulevement des particules dans le modele convectif
2378!c
2379      do i = 1,klon
2380        ALE(i) = 0.
2381        ALP(i) = 0.
2382      enddo
2383!c
2384!ccalcul de ale_wake et alp_wake
2385       if (iflag_wake>=1) then
2386         if (itap .le. it_wape_prescr) then
2387          do i = 1,klon
2388           ale_wake(i) = wape_prescr
2389           alp_wake(i) = fip_prescr
2390          enddo
2391         else
2392          do i = 1,klon
2393!cjyg  ALE=WAPE au lieu de ALE = 1/2 Cstar**2
2394!ccc           ale_wake(i) = 0.5*wake_cstar(i)**2
2395           ale_wake(i) = wake_pe(i)
2396           alp_wake(i) = wake_fip(i)
2397          enddo
2398         endif
2399       else
2400         do i = 1,klon
2401           ale_wake(i) = 0.
2402           alp_wake(i) = 0.
2403         enddo
2404       endif
2405!ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2406!cdans le thermique sinon
2407       if (iflag_coupl.eq.0) then
2408          if (debut.and.prt_level.gt.9)                                              &
2409       &                     WRITE(lunout,*)'ALE et ALP imposes'
2410          do i = 1,klon
2411!con ne couple que ale
2412!c           ALE(i) = max(ale_wake(i),Ale_bl(i))
2413            ALE(i) = max(ale_wake(i),ale_bl_prescr)
2414!con ne couple que alp
2415!c           ALP(i) = alp_wake(i) + Alp_bl(i)
2416            ALP(i) = alp_wake(i) + alp_bl_prescr
2417          enddo
2418       else
2419         IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2420!         do i = 1,klon
2421!             ALE(i) = max(ale_wake(i),Ale_bl(i))
2422! avant        ALP(i) = alp_wake(i) + Alp_bl(i)
2423!             ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2424!         write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2425!         write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2426!         enddo
2427
2428!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2429! Modif FH 2010/04/27. Sans doute temporaire.
2430! Deux options pour le alp_offset : constant si >?? 0 ou proportionnel ??a
2431! w si <0
2432!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2433       do i = 1,klon
2434          ALE(i) = max(ale_wake(i),Ale_bl(i))
2435!ccc nrlmd le 10/04/2012----------Stochastic triggering--------------
2436          if (iflag_trig_bl.ge.1) then
2437             ALE(i) = max(ale_wake(i),Ale_bl_trig(i))
2438          endif
2439!ccc fin nrlmd le 10/04/2012
2440          if (alp_offset>=0.) then
2441            ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2442          else
2443            ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
2444            if (alp(i)<0.) then
2445               print*,'ALP ',alp(i),alp_wake(i)                                      &
2446       &         ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
2447            endif
2448          endif
2449       enddo
2450!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2451
2452       endif
2453       do i=1,klon
2454          if (alp(i)>alp_max) then
2455             IF(prt_level>9)WRITE(lunout,*)                             &
2456       &       'WARNING SUPER ALP (seuil=',alp_max,                                  &
2457       &       '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2458             alp(i)=alp_max
2459          endif
2460          if (ale(i)>ale_max) then
2461             IF(prt_level>9)WRITE(lunout,*)                             &
2462       &       'WARNING SUPER ALE (seuil=',ale_max,                                  &
2463       &       '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2464             ale(i)=ale_max
2465          endif
2466       enddo
2467
2468!cfin calcul ale et alp
2469!c=================================================================================================
2470
2471!c sb, oct02:
2472!c Schema de convection modularise et vectorise:
2473!c (driver commun aux versions 3 et 4)
2474!c
2475          IF (ok_cvl) THEN ! new driver for convectL
2476
2477             IF (type_trac == 'repr') THEN
2478                nbtr_tmp=ntra
2479             ELSE
2480                nbtr_tmp=nbtr
2481             END IF
2482          CALL concvl (iflag_con,iflag_clos,                                         &
2483       &        dtime,paprs,pplay,t_undi,q_undi,                                     &
2484       &        t_wake,q_wake,wake_s,                                                &
2485       &        u_seri,v_seri,tr_seri,nbtr_tmp,                                      &
2486       &        ALE,ALP,                                                             &
2487       &        ema_work1,ema_work2,                                                 &
2488       &        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,                                &
2489       &        rain_con, snow_con, ibas_con, itop_con, sigd,                        &
2490       &        ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0,                            &
2491       &        Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,                         &
2492       &        pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,               &
2493! RomP >>>
2494!!     .        pmflxr,pmflxs,da,phi,mp,
2495!!     .        ftd,fqd,lalim_conv,wght_th)                                          &
2496       &        pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij,                   &
2497       &        ftd,fqd,lalim_conv,wght_th,                                          &
2498       &        ev, ep,epmlmMm,eplaMm,                                               &
2499       &        wdtrainA,wdtrainM)
2500! RomP <<<
2501
2502!IM begin
2503!c       print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2504!c    .dnwd0(1,1),ftd(1,1),fqd(1,1)
2505!IM end
2506!IM cf. FH
2507              clwcon0=qcondc
2508              pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2509
2510              do i = 1, klon
2511                if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2512              enddo
2513
2514          ELSE ! ok_cvl
2515
2516!c MAF conema3 ne contient pas les traceurs
2517          CALL conema3 (dtime,                                                       &
2518       &        paprs,pplay,t_seri,q_seri,                                           &
2519       &        u_seri,v_seri,tr_seri,ntra,                                          &
2520       &        ema_work1,ema_work2,                                                 &
2521       &        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,                                &
2522       &        rain_con, snow_con, ibas_con, itop_con,                              &
2523       &        upwd,dnwd,dnwd0,bas,top,                                             &
2524       &        Ma,cape,tvp,rflag,                                                   &
2525       &        pbase                                                                &
2526       &        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr                               &
2527       &        ,clwcon0)
2528
2529          ENDIF ! ok_cvl
2530
2531!c
2532!c Correction precip
2533          rain_con = rain_con * cvl_corr
2534          snow_con = snow_con * cvl_corr
2535!c
2536
2537           IF (.NOT. ok_gust) THEN
2538           do i = 1, klon
2539            wd(i)=0.0
2540           enddo
2541           ENDIF
2542
2543!c =================================================================== c
2544!c Calcul des proprietes des nuages convectifs
2545!c
2546
2547!c   calcul des proprietes des nuages convectifs
2548             clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2549             call clouds_gno                                                         &
2550       &       (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2551
2552!c =================================================================== c
2553
2554          DO i = 1, klon
2555           itop_con(i) = min(max(itop_con(i),1),klev)
2556           ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2557          ENDDO
2558
2559          DO i = 1, klon
2560            ema_pcb(i)  = paprs(i,ibas_con(i))
2561          ENDDO
2562          DO i = 1, klon
2563! L'idicage de itop_con peut cacher un pb potentiel
2564! FH sous la dictee de JYG, CR
2565            ema_pct(i)  = paprs(i,itop_con(i)+1)
2566
2567            if (itop_con(i).gt.klev-3) then
2568              if(prt_level >= 9) then
2569                write(lunout,*)'La convection monte trop haut '
2570                write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2571              endif
2572            endif
2573          ENDDO     
2574      ELSE IF (iflag_con.eq.0) THEN
2575          write(lunout,*) 'On n appelle pas la convection'
2576          clwcon0=0.
2577          rnebcon0=0.
2578          d_t_con=0.
2579          d_q_con=0.
2580          d_u_con=0.
2581          d_v_con=0.
2582          rain_con=0.
2583          snow_con=0.
2584          bas=1
2585          top=1
2586      ELSE
2587          WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2588          CALL abort
2589      ENDIF
2590
2591!c     CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2592!c    .              d_u_con, d_v_con)
2593
2594!-----------------------------------------------------------------------------------------
2595! ajout des tendances de la diffusion turbulente
2596      CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2597!-----------------------------------------------------------------------------------------
2598
2599      if (mydebug) then
2600        call writefield_phy('u_seri',u_seri,llm)
2601        call writefield_phy('v_seri',v_seri,llm)
2602        call writefield_phy('t_seri',t_seri,llm)
2603        call writefield_phy('q_seri',q_seri,llm)
2604      endif
2605
2606!IM
2607      IF (ip_ebil_phy.ge.2) THEN
2608        ztit='after convect'
2609        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime                             &
2610       &      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay              &
2611       &      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2612         call diagphy(airephy,ztit,ip_ebil_phy                                       &
2613       &      , zero_v, zero_v, zero_v, zero_v, zero_v                               &
2614       &      , zero_v, rain_con, snow_con, ztsol                                    &
2615       &      , d_h_vcol, d_qt, d_ec                                                 &
2616       &      , fs_bound, fq_bound )
2617      END IF
2618!C
2619      PRINT *,'  Lluis 4 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
2620      IF (check) THEN
2621          za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2622          WRITE(lunout,*)"aprescon=", za
2623          zx_t = 0.0
2624          za = 0.0
2625          DO i = 1, klon
2626            za = za + airephy(i)/REAL(klon)
2627            zx_t = zx_t + (rain_con(i)+                                              &
2628       &                   snow_con(i))*airephy(i)/REAL(klon)
2629          ENDDO
2630          zx_t = zx_t/za*dtime
2631          WRITE(lunout,*)"Precip=", zx_t
2632      ENDIF
2633      IF (zx_ajustq) THEN
2634          DO i = 1, klon
2635            z_apres(i) = 0.0
2636          ENDDO
2637          DO k = 1, klev
2638            DO i = 1, klon
2639              z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))                   &
2640       &            *(paprs(i,k)-paprs(i,k+1))/RG
2641            ENDDO
2642          ENDDO
2643          DO i = 1, klon
2644            z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)               &
2645       &          /z_apres(i)
2646          ENDDO
2647          DO k = 1, klev
2648            DO i = 1, klon
2649              IF (z_factor(i).GT.(1.0+1.0E-08) .OR.                                  &
2650       &            z_factor(i).LT.(1.0-1.0E-08)) THEN
2651                  q_seri(i,k) = q_seri(i,k) * z_factor(i)
2652              ENDIF
2653            ENDDO
2654          ENDDO
2655      ENDIF
2656      zx_ajustq=.FALSE.
2657
2658!c
2659!c=============================================================================
2660!cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2661!cpour la couche limite diffuse pour l instant
2662!c
2663      if (iflag_wake>=1) then
2664      DO k=1,klev
2665        DO i=1,klon
2666          dt_dwn(i,k)  = ftd(i,k)
2667          wdt_PBL(i,k) = 0.
2668          dq_dwn(i,k)  = fqd(i,k)
2669          wdq_PBL(i,k) = 0.
2670          M_dwn(i,k)   = dnwd0(i,k)
2671          M_up(i,k)    = upwd(i,k)
2672          dt_a(i,k)    = d_t_con(i,k)/dtime - ftd(i,k)
2673          udt_PBL(i,k) = 0.
2674          dq_a(i,k)    = d_q_con(i,k)/dtime - fqd(i,k)
2675          udq_PBL(i,k) = 0.
2676        ENDDO
2677      ENDDO
2678
2679      if (iflag_wake==2) then
2680        ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2681        DO k = 1,klev
2682         dt_dwn(:,k)= dt_dwn(:,k)+                                                   &
2683       &            ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
2684         dq_dwn(:,k)= dq_dwn(:,k)+                                                   &
2685       &            ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
2686        ENDDO
2687      endif
2688!c
2689!ccalcul caracteristiques de la poche froide
2690      call calWAKE (paprs,pplay,dtime                                                &
2691       &               ,t_seri,q_seri,omega                                          &
2692       &               ,dt_dwn,dq_dwn,M_dwn,M_up                                     &
2693       &               ,dt_a,dq_a,sigd                                               &
2694       &               ,wdt_PBL,wdq_PBL                                              &
2695       &               ,udt_PBL,udq_PBL                                              &
2696       &               ,wake_deltat,wake_deltaq,wake_dth                             &
2697       &               ,wake_h,wake_s,wake_dens                                      &
2698       &               ,wake_pe,wake_fip,wake_gfl                                    &
2699       &               ,dt_wake,dq_wake                                              &
2700       &               ,wake_k, t_undi,q_undi                                        &
2701       &               ,wake_omgbdth,wake_dp_omgb                                    &
2702       &               ,wake_dtKE,wake_dqKE                                          &
2703       &               ,wake_dtPBL,wake_dqPBL                                        &
2704       &               ,wake_omg,wake_dp_deltomg                                     &
2705       &               ,wake_spread,wake_Cstar,wake_d_deltat_gw                      &
2706       &               ,wake_ddeltat,wake_ddeltaq)
2707!c
2708!-----------------------------------------------------------------------------------------
2709! ajout des tendances des poches froides
2710! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2711! coherent avec les autres d_t_...
2712      d_t_wake(:,:)=dt_wake(:,:)*dtime
2713      d_q_wake(:,:)=dq_wake(:,:)*dtime
2714      CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2715!-----------------------------------------------------------------------------------------
2716
2717      endif
2718!c
2719!c===================================================================
2720!cJYG
2721      IF (ip_ebil_phy.ge.2) THEN
2722        ztit='after wake'
2723        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime                             &
2724       &      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay              &
2725       &      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2726        call diagphy(airephy,ztit,ip_ebil_phy                                        &
2727       &      , zero_v, zero_v, zero_v, zero_v, zero_v                               &
2728       &      , zero_v, zero_v, zero_v, ztsol                                        &
2729       &      , d_h_vcol, d_qt, d_ec                                                 &
2730       &      , fs_bound, fq_bound )
2731      END IF
2732      PRINT *,'  Lluis 5 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
2733!c      print*,'apres callwake iflag_cldcon=', iflag_cldcon
2734!c
2735!c===================================================================
2736!c Convection seche (thermiques ou ajustement)
2737!c===================================================================
2738!c
2739       call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri                         &
2740       & ,seuil_inversion,weak_inversion,dthmin)
2741
2742
2743
2744      d_t_ajsb(:,:)=0.
2745      d_q_ajsb(:,:)=0.
2746      d_t_ajs(:,:)=0.
2747      d_u_ajs(:,:)=0.
2748      d_v_ajs(:,:)=0.
2749      d_q_ajs(:,:)=0.
2750      clwcon0th(:,:)=0.
2751!c
2752!c      fm_therm(:,:)=0.
2753!c      entr_therm(:,:)=0.
2754!c      detr_therm(:,:)=0.
2755!c
2756
2757      IF(prt_level>9)WRITE(lunout,*)                                                 &
2758       &    'AVANT LA CONVECTION SECHE , iflag_thermals='                            &
2759       &   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
2760      if(iflag_thermals.lt.0) then
2761!c  Rien
2762!c  ====
2763         IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2764
2765
2766      else
2767
2768!c  Thermiques
2769!c  ==========
2770         IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='               &
2771       &   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
2772
2773
2774!ccc nrlmd le 10/04/2012
2775         DO k=1,klev+1
2776           DO i=1,klon
2777           pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
2778           pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
2779           pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
2780           pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
2781           ENDDO
2782         ENDDO
2783!ccc fin nrlmd le 10/04/2012
2784
2785         if (iflag_thermals>=1) then
2786         call calltherm(pdtphys                                                      &
2787       &      ,pplay,paprs,pphi,weak_inversion                                       &
2788       &      ,u_seri,v_seri,t_seri,q_seri,zqsat,debut                               &
2789       &      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs                                       &
2790       &      ,fm_therm,entr_therm,detr_therm                                        &
2791       &      ,zqasc,clwcon0th,lmax_th,ratqscth                                      &
2792       &      ,ratqsdiff,zqsatth                                                     &
2793!con rajoute ale et alp, et les caracteristiques de la couche alim                   &
2794       &      ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca                &
2795       &      ,ztv,zpspsk,ztla,zthl                                                  &
2796!ccc nrlmd le 10/04/2012                                                             &
2797       &      ,pbl_tke_input,pctsrf,omega,airephy                                    &
2798       &      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0                  &
2799       &      ,n2,s2,ale_bl_stat                                                     &
2800       &      ,therm_tke_max,env_tke_max                                             &
2801       &      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke                            &
2802       &      ,alp_bl_conv,alp_bl_stat                                               &
2803!ccc fin nrlmd le 10/04/2012                                                         &
2804       &      ,zqla,ztva )
2805
2806!ccc nrlmd le 10/04/2012
2807!c-----------Stochastic triggering-----------
2808      if (iflag_trig_bl.ge.1) then
2809!c
2810        IF (prt_level .GE. 10) THEN
2811         print *,'cin, ale_bl_stat, alp_bl_stat ',                                   &
2812       &            cin, ale_bl_stat, alp_bl_stat
2813        ENDIF
2814
2815!c----Initialisations
2816      do i=1,klon
2817      proba_notrig(i)=1.
2818      random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
2819        if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
2820        tau_trig(i)=tau_trig_shallow
2821        else
2822        tau_trig(i)=tau_trig_deep
2823        endif
2824      enddo
2825!c
2826        IF (prt_level .GE. 10) THEN
2827         print *,'random_notrig, tau_trig ',                                         &
2828       &            random_notrig, tau_trig
2829          print *,'s_trig,s2,n2 ',                                                   &
2830       &             s_trig,s2,n2
2831        ENDIF
2832
2833!c----Tirage al\'eatoire et calcul de ale_bl_trig
2834      do i=1,klon
2835        if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) )  then
2836        proba_notrig(i)=(1.-exp(-s_trig/s2(i)))**                                    &
2837       &                  (n2(i)*dtime/tau_trig(i))
2838!c        print *, 'proba_notrig(i) ',proba_notrig(i)
2839          if (random_notrig(i) .ge. proba_notrig(i)) then
2840          ale_bl_trig(i)=ale_bl_stat(i)
2841          else
2842          ale_bl_trig(i)=0.
2843          endif
2844        else
2845        proba_notrig(i)=1.
2846        random_notrig(i)=0.
2847        ale_bl_trig(i)=0.
2848        endif
2849      enddo
2850!c
2851        IF (prt_level .GE. 10) THEN
2852         print *,'proba_notrig, ale_bl_trig ',                                       &
2853       &            proba_notrig, ale_bl_trig
2854        ENDIF
2855
2856      endif !(iflag_trig_bl)
2857
2858!c-----------Statistical closure-----------
2859      if (iflag_clos_bl.ge.1) then
2860
2861        do i=1,klon
2862        alp_bl(i)=alp_bl_stat(i)
2863        enddo
2864
2865      else
2866
2867      alp_bl_stat(:)=0.
2868
2869      endif !(iflag_clos_bl)
2870
2871        IF (prt_level .GE. 10) THEN
2872         print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
2873        ENDIF
2874
2875!ccc fin nrlmd le 10/04/2012
2876
2877! ----------------------------------------------------------------------
2878! Transport de la TKE par les panaches thermiques.
2879! FH : 2010/02/01
2880!     if (iflag_pbl.eq.10) then
2881!     call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
2882!    s           rg,paprs,pbl_tke)
2883!     endif
2884! ----------------------------------------------------------------------
2885!IM/FH: 2011/02/23
2886! Couplage Thermiques/Emanuel seulement si T<0
2887      if (iflag_coupl==2) then
2888        print*,'Couplage Thermiques/Emanuel seulement si T<0'
2889        do i=1,klon
2890           if (t_seri(i,lmax_th(i))>273.) then
2891              Ale_bl(i)=0.
2892           endif
2893        enddo
2894      endif
2895
2896      do i=1,klon
2897         zmax_th(i)=pphi(i,lmax_th(i))/rg
2898      enddo
2899
2900         endif
2901
2902
2903!c  Ajustement sec
2904!c  ==============
2905
2906! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
2907! a partir du sommet des thermiques.
2908! Dans le cas contraire, on demarre au niveau 1.
2909
2910         if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2911
2912         if(iflag_thermals.eq.0) then
2913            IF(prt_level>9)WRITE(lunout,*)'ajsec'
2914            limbas(:)=1
2915         else
2916            limbas(:)=lmax_th(:)
2917         endif
2918 
2919! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2920! pour des test de convergence numerique.
2921! Le nouveau ajsec est a priori mieux, meme pour le cas
2922! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2923! non nulles numeriquement pour des mailles non concernees.
2924
2925         if (iflag_thermals.eq.0) then
2926            CALL ajsec_convV2(paprs, pplay, t_seri,q_seri                            &
2927       &      , d_t_ajsb, d_q_ajsb)
2928         else
2929            CALL ajsec(paprs, pplay, t_seri,q_seri,limbas                            &
2930       &      , d_t_ajsb, d_q_ajsb)
2931         endif
2932
2933!-----------------------------------------------------------------------------------------
2934! ajout des tendances de l'ajustement sec ou des thermiques
2935      CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2936         d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2937         d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2938
2939!-----------------------------------------------------------------------------------------
2940
2941         endif
2942
2943      endif
2944
2945!c
2946!c===================================================================
2947!IM
2948      IF (ip_ebil_phy.ge.2) THEN
2949        ztit='after dry_adjust'
2950        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime                             &
2951       &      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay              &
2952       &      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2953        call diagphy(airephy,ztit,ip_ebil_phy                                        &
2954       &      , zero_v, zero_v, zero_v, zero_v, zero_v                               &
2955       &      , zero_v, zero_v, zero_v, ztsol                                        &
2956       &      , d_h_vcol, d_qt, d_ec                                                 &
2957       &      , fs_bound, fq_bound )
2958      END IF
2959      PRINT *,'  Lluis 6 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
2960
2961!c-------------------------------------------------------------------------
2962! Computation of ratqs, the width (normalized) of the subrid scale
2963! water distribution
2964      CALL  calcratqs(klon,klev,prt_level,lunout,                                    &
2965       &     iflag_ratqs,iflag_con,iflag_cldcon,pdtphys,                             &
2966       &     ratqsbas,ratqshaut,tau_ratqs,fact_cldcon,                               &
2967       &     ptconv,ptconvth,clwcon0th, rnebcon0th,                                  &
2968       &     paprs,pplay,q_seri,zqsat,fm_therm,                                      &
2969       &     ratqs,ratqsc)
2970
2971!c
2972!c Appeler le processus de condensation a grande echelle
2973!c et le processus de precipitation
2974!c-------------------------------------------------------------------------
2975      IF (prt_level .GE.10) THEN
2976       print *,' ->fisrtilp '
2977      ENDIF
2978!c-------------------------------------------------------------------------
2979      CALL fisrtilp(dtime,paprs,pplay,                                               &
2980       &           t_seri, q_seri,ptconv,ratqs,                                      &
2981       &           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,                         &
2982       &           rain_lsc, snow_lsc,                                               &
2983       &           pfrac_impa, pfrac_nucl, pfrac_1nucl,                              &
2984       &           frac_impa, frac_nucl, beta_prec_fisrt,                            &
2985       &           prfl, psfl, rhcl,                                                 &
2986       &           zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon )
2987
2988      WHERE (rain_lsc < 0) rain_lsc = 0.
2989      WHERE (snow_lsc < 0) snow_lsc = 0.
2990!-----------------------------------------------------------------------------------------
2991! ajout des tendances de la diffusion turbulente
2992      CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
2993!-----------------------------------------------------------------------------------------
2994      DO k = 1, klev
2995      DO i = 1, klon
2996         cldfra(i,k) = rneb(i,k)
2997         IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2998      ENDDO
2999      ENDDO
3000      IF (check) THEN
3001         za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
3002         WRITE(lunout,*)"apresilp=", za
3003         zx_t = 0.0
3004         za = 0.0
3005         DO i = 1, klon
3006            za = za + airephy(i)/REAL(klon)
3007            zx_t = zx_t + (rain_lsc(i)                                               &
3008       &                  + snow_lsc(i))*airephy(i)/REAL(klon)
3009        ENDDO
3010         zx_t = zx_t/za*dtime
3011         WRITE(lunout,*)"Precip=", zx_t
3012      ENDIF
3013!IM
3014      IF (ip_ebil_phy.ge.2) THEN
3015        ztit='after fisrt'
3016        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime                             &
3017       &      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay              &
3018       &      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3019        call diagphy(airephy,ztit,ip_ebil_phy                                        &
3020       &      , zero_v, zero_v, zero_v, zero_v, zero_v                               &
3021       &      , zero_v, rain_lsc, snow_lsc, ztsol                                    &
3022       &      , d_h_vcol, d_qt, d_ec                                                 &
3023       &      , fs_bound, fq_bound )
3024      END IF
3025      PRINT *,'  Lluis 7 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
3026      if (mydebug) then
3027        call writefield_phy('u_seri',u_seri,llm)
3028        call writefield_phy('v_seri',v_seri,llm)
3029        call writefield_phy('t_seri',t_seri,llm)
3030        call writefield_phy('q_seri',q_seri,llm)
3031      endif
3032!c
3033!c-------------------------------------------------------------------
3034!c  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3035!c-------------------------------------------------------------------
3036
3037!c 1. NUAGES CONVECTIFS
3038!c
3039!IM cf FH
3040!c     IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
3041      IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
3042       snow_tiedtke=0.
3043!c     print*,'avant calcul de la pseudo precip '
3044!c     print*,'iflag_cldcon',iflag_cldcon
3045       if (iflag_cldcon.eq.-1) then
3046          rain_tiedtke=rain_con
3047       else
3048!c       print*,'calcul de la pseudo precip '
3049          rain_tiedtke=0.
3050!c         print*,'calcul de la pseudo precip 0'
3051          do k=1,klev
3052          do i=1,klon
3053             if (d_q_con(i,k).lt.0.) then
3054                rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys                 &
3055       &         *(paprs(i,k)-paprs(i,k+1))/rg
3056             endif
3057          enddo
3058          enddo
3059       endif
3060!c
3061!c     call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3062!c
3063
3064!c Nuages diagnostiques pour Tiedtke
3065      CALL diagcld1(paprs,pplay,                                                     &
3066!IM cf FH  .             rain_con,snow_con,ibas_con,itop_con,                       &
3067       &             rain_tiedtke,snow_tiedtke,ibas_con,itop_con,                    &
3068       &             diafra,dialiq)
3069      DO k = 1, klev
3070      DO i = 1, klon
3071      IF (diafra(i,k).GT.cldfra(i,k)) THEN
3072         cldliq(i,k) = dialiq(i,k)
3073         cldfra(i,k) = diafra(i,k)
3074      ENDIF
3075      ENDDO
3076      ENDDO
3077
3078      ELSE IF (iflag_cldcon.ge.3) THEN
3079!c  On prend pour les nuages convectifs le max du calcul de la
3080!c  convection et du calcul du pas de temps precedent diminue d'un facteur
3081!c  facttemps
3082      facteur = pdtphys *facttemps
3083      do k=1,klev
3084         do i=1,klon
3085            rnebcon(i,k)=rnebcon(i,k)*facteur
3086            if (rnebcon0(i,k)*clwcon0(i,k).gt.rnebcon(i,k)*clwcon(i,k))              &
3087       &      then
3088                rnebcon(i,k)=rnebcon0(i,k)
3089                clwcon(i,k)=clwcon0(i,k)
3090            endif
3091         enddo
3092      enddo
3093
3094!c
3095!cjq - introduce the aerosol direct and first indirect radiative forcings
3096!cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
3097      IF (flag_aerosol .gt. 0) THEN
3098         IF (.NOT. aerosol_couple)                                                   &
3099       &        CALL readaerosol_optic(                                              &
3100       &        debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref,                   &
3101       &        pdtphys, pplay, paprs, t_seri, rhcl, presnivs,                       &
3102       &        mass_solu_aero, mass_solu_aero_pi,                                   &
3103       &        tau_aero, piz_aero, cg_aero,                                         &
3104       &        tausum_aero, tau3d_aero)
3105      ELSE
3106         tausum_aero(:,:,:) = 0.
3107         tau_aero(:,:,:,:) = 0.
3108         piz_aero(:,:,:,:) = 0.
3109         cg_aero(:,:,:,:)  = 0.
3110      ENDIF
3111!c
3112!c--STRAT AEROSOL
3113!c--updates tausum_aero,tau_aero,piz_aero,cg_aero
3114      IF (flag_aerosol_strat) THEN
3115         PRINT *,'appel a readaerosolstrat', mth_cur
3116         CALL readaerosolstrato(debut)
3117      ENDIF
3118!c--fin STRAT AEROSOL
3119
3120!IM calcul nuages par le simulateur ISCCP
3121!c
3122#ifdef histISCCP
3123      IF (ok_isccp) THEN
3124!c
3125!IM lecture invtau, tautab des fichiers formattes
3126!c
3127      IF (debut) THEN
3128!$OMP MASTER
3129!c
3130      open(99,file='tautab.formatted', FORM='FORMATTED')
3131      read(99,'(f30.20)') tautab_omp
3132      close(99)
3133!c
3134      open(99,file='invtau.formatted',form='FORMATTED')
3135      read(99,'(i10)') invtau_omp
3136
3137!c     print*,'calcul_simulISCCP invtau_omp',invtau_omp
3138!c     write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
3139
3140      close(99)
3141!$OMP END MASTER
3142!$OMP BARRIER
3143      tautab=tautab_omp
3144      invtau=invtau_omp
3145!c
3146      ENDIF !debut
3147!c
3148!IM appel simulateur toutes les  NINT(freq_ISCCP/dtime) heures
3149       IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
3150#include "calcul_simulISCCP.h"
3151       ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
3152      ENDIF !ok_isccp
3153#endif
3154
3155!c   On prend la somme des fractions nuageuses et des contenus en eau
3156
3157      if (iflag_cldcon>=5) then
3158
3159        do k=1,klev
3160         ptconvth(:,k)=fm_therm(:,k+1)>0.
3161        enddo
3162
3163       if (iflag_coupl==4) then
3164
3165! Dans le cas iflag_coupl==4, on prend la somme des convertures
3166! convectives et lsc dans la partie des thermiques
3167! Le controle par iflag_coupl est peut etre provisoire.
3168         do k=1,klev
3169            do i=1,klon
3170               if (ptconv(i,k).and.ptconvth(i,k)) then
3171                   cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3172                   cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3173               else if (ptconv(i,k)) then
3174                   cldfra(i,k)=rnebcon(i,k)
3175                   cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3176               endif
3177            enddo
3178         enddo
3179
3180         else if (iflag_coupl==5) then
3181         do k=1,klev
3182            do i=1,klon
3183               cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3184               cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3185            enddo
3186         enddo
3187
3188         else
3189
3190! Si on est sur un point touche par la convection profonde et pas
3191! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
3192! de la convection profonde.
3193
3194!IM/FH: 2011/02/23
3195! definition des points sur lesquels ls thermiques sont actifs
3196
3197         do k=1,klev
3198            do i=1,klon
3199               if (ptconv(i,k).and. .not. ptconvth(i,k)) then
3200                   cldfra(i,k)=rnebcon(i,k)
3201                   cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3202               endif
3203            enddo
3204         enddo
3205
3206        endif
3207
3208      else
3209
3210! Ancienne version
3211      cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3212      cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3213      endif
3214
3215      ENDIF
3216
3217!     plulsc(:)=0.
3218!     do k=1,klev,-1
3219!        do i=1,klon
3220!              zzz=prfl(:,k)+psfl(:,k)
3221!           if (.not.ptconvth.zzz.gt.0.)
3222!        enddo prfl, psfl,
3223!     enddo
3224!c
3225!c 2. NUAGES STARTIFORMES
3226!c
3227      IF (ok_stratus) THEN
3228      CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3229      DO k = 1, klev
3230      DO i = 1, klon
3231      IF (diafra(i,k).GT.cldfra(i,k)) THEN
3232         cldliq(i,k) = dialiq(i,k)
3233         cldfra(i,k) = diafra(i,k)
3234      ENDIF
3235      ENDDO
3236      ENDDO
3237      ENDIF
3238!c
3239!c Precipitation totale
3240!c
3241      DO i = 1, klon
3242         rain_fall(i) = rain_con(i) + rain_lsc(i)
3243         snow_fall(i) = snow_con(i) + snow_lsc(i)
3244      ENDDO
3245!IM
3246      IF (ip_ebil_phy.ge.2) THEN
3247        ztit="after diagcld"
3248        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime                             &
3249       &      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay              &
3250       &      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3251        call diagphy(airephy,ztit,ip_ebil_phy                                        &
3252       &      , zero_v, zero_v, zero_v, zero_v, zero_v                               &
3253       &      , zero_v, zero_v, zero_v, ztsol                                        &
3254       &      , d_h_vcol, d_qt, d_ec                                                 &
3255       &      , fs_bound, fq_bound )
3256      END IF
3257      PRINT *,'  Lluis 8 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
3258!c
3259!c Calculer l'humidite relative pour diagnostique
3260!c
3261      DO k = 1, klev
3262      DO i = 1, klon
3263         zx_t = t_seri(i,k)
3264         IF (thermcep) THEN
3265            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3266            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3267            zx_qs  = MIN(0.5,zx_qs)
3268            zcor   = 1./(1.-retv*zx_qs)
3269            zx_qs  = zx_qs*zcor
3270         ELSE
3271           IF (zx_t.LT.t_coup) THEN
3272              zx_qs = qsats(zx_t)/pplay(i,k)
3273           ELSE
3274              zx_qs = qsatl(zx_t)/pplay(i,k)
3275           ENDIF
3276         ENDIF
3277         zx_rh(i,k) = q_seri(i,k)/zx_qs
3278         zqsat(i,k)=zx_qs
3279      ENDDO
3280      ENDDO
3281
3282!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3283!c   equivalente a 2m (tpote) pour diagnostique
3284!c
3285      DO i = 1, klon
3286       tpot(i)=zt2m(i)*(100000./paprs(i,1))**RKAPPA
3287       IF (thermcep) THEN
3288        IF(zt2m(i).LT.RTT) then
3289         Lheat=RLSTT
3290        ELSE
3291         Lheat=RLVTT
3292        ENDIF
3293       ELSE
3294        IF (zt2m(i).LT.RTT) THEN
3295         Lheat=RLSTT
3296        ELSE
3297         Lheat=RLVTT
3298        ENDIF
3299       ENDIF
3300       tpote(i) = tpot(i)*                                                           &
3301       & EXP((Lheat *qsat2m(i))/(RCPD*zt2m(i)))
3302      ENDDO
3303
3304      IF (type_trac == 'inca') THEN
3305#ifdef INCA
3306         CALL VTe(VTphysiq)
3307         CALL VTb(VTinca)
3308         calday = REAL(days_elapsed + 1) + jH_cur
3309
3310         call chemtime(itap+itau_phy-1, date0, dtime)
3311         IF (config_inca == 'aero') THEN
3312            CALL AEROSOL_METEO_CALC(                                                 &
3313       &           calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,                       &
3314       &           prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
3315         END IF
3316
3317         zxsnow_dummy(:) = 0.0
3318
3319         CALL chemhook_begin (calday,                                                &
3320       &                          days_elapsed+1,                                    &
3321       &                          jH_cur,                                            &
3322       &                          pctsrf(1,1),                                       &
3323       &                          rlat,                                              &
3324       &                          rlon,                                              &
3325       &                          airephy,                                           &
3326       &                          paprs,                                             &
3327       &                          pplay,                                             &
3328       &                          coefh(:,:,is_ave),                                 &
3329       &                          pphi,                                              &
3330       &                          t_seri,                                            &
3331       &                          u,                                                 &
3332       &                          v,                                                 &
3333       &                          wo(:, :, 1),                                       &
3334       &                          q_seri,                                            &
3335       &                          zxtsol,                                            &
3336       &                          zxsnow_dummy,                                      &
3337       &                          solsw,                                             &
3338       &                          albsol1,                                           &
3339       &                          rain_fall,                                         &
3340       &                          snow_fall,                                         &
3341       &                          itop_con,                                          &
3342       &                          ibas_con,                                          &
3343       &                          cldfra,                                            &
3344       &                          iim,                                               &
3345       &                          jjm,                                               &
3346       &                          tr_seri,                                           &
3347       &                          ftsol,                                             &
3348       &                          paprs,                                             &
3349       &                          cdragh,                                            &
3350       &                          cdragm,                                            &
3351       &                          pctsrf,                                            &
3352       &                          pdtphys,                                           &
3353       &                            itap)
3354
3355         CALL VTe(VTinca)
3356         CALL VTb(VTphysiq)
3357#endif
3358      END IF !type_trac = inca
3359!c     
3360!c Calculer les parametres optiques des nuages et quelques
3361!c parametres pour diagnostiques:
3362!c
3363
3364      IF (aerosol_couple) THEN
3365         mass_solu_aero(:,:)    = ccm(:,:,1)
3366         mass_solu_aero_pi(:,:) = ccm(:,:,2)
3367      END IF
3368
3369      if (ok_newmicro) then
3370      CALL newmicro (ok_cdnc, bl95_b0, bl95_b1,                                      &
3371       &              paprs, pplay, t_seri, cldliq, cldfra,                          &
3372       &              cldtau, cldemi, cldh, cldl, cldm, cldt, cldq,                  &
3373       &              flwp, fiwp, flwc, fiwc,                                        &
3374       &              mass_solu_aero, mass_solu_aero_pi,                             &
3375       &              cldtaupi, re, fl, ref_liq, ref_ice)
3376      else
3377      CALL nuage (paprs, pplay,                                                      &
3378       &            t_seri, cldliq, cldfra, cldtau, cldemi,                          &
3379       &            cldh, cldl, cldm, cldt, cldq,                                    &
3380       &            ok_aie,                                                          &
3381       &            mass_solu_aero, mass_solu_aero_pi,                               &
3382       &            bl95_b0, bl95_b1,                                                &
3383       &            cldtaupi, re, fl)
3384      endif
3385!c
3386!IM betaCRF
3387!c
3388      cldtaurad   = cldtau
3389      cldtaupirad = cldtaupi
3390      cldemirad   = cldemi                                                           &
3391       &!c
3392      if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND.                               &
3393       &lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3394!c
3395!c global
3396!c
3397       DO k=1, klev
3398       DO i=1, klon
3399        if (pplay(i,k).GE.pfree) THEN
3400         beta(i,k) = beta_pbl
3401        else
3402         beta(i,k) = beta_free
3403        endif
3404        if (mskocean_beta) THEN
3405         beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3406        endif
3407        cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
3408        cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3409        cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
3410        cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
3411       ENDDO
3412       ENDDO
3413!c
3414      else
3415!c
3416!c regional
3417!c
3418       DO k=1, klev
3419       DO i=1,klon
3420!c
3421        if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND.                       &
3422       &      rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
3423         if (pplay(i,k).GE.pfree) THEN
3424          beta(i,k) = beta_pbl
3425         else
3426          beta(i,k) = beta_free
3427         endif
3428         if (mskocean_beta) THEN
3429          beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3430         endif
3431        cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
3432        cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3433        cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
3434        cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
3435        endif
3436!c
3437       ENDDO
3438       ENDDO
3439!c
3440      endif
3441!c
3442!c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3443!c
3444      IF (MOD(itaprad,radpas).EQ.0) THEN
3445
3446      DO i = 1, klon
3447         albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)                             &
3448       &             + falb1(i,is_lic) * pctsrf(i,is_lic)                            &
3449       &             + falb1(i,is_ter) * pctsrf(i,is_ter)                            &
3450       &             + falb1(i,is_sic) * pctsrf(i,is_sic)
3451         albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)                             &
3452       &               + falb2(i,is_lic) * pctsrf(i,is_lic)                          &
3453       &               + falb2(i,is_ter) * pctsrf(i,is_ter)                          &
3454       &               + falb2(i,is_sic) * pctsrf(i,is_sic)
3455      ENDDO
3456
3457      if (mydebug) then
3458        call writefield_phy('u_seri',u_seri,llm)
3459        call writefield_phy('v_seri',v_seri,llm)
3460        call writefield_phy('t_seri',t_seri,llm)
3461        call writefield_phy('q_seri',q_seri,llm)
3462      endif
3463     
3464      IF (aerosol_couple) THEN
3465#ifdef INCA
3466         CALL radlwsw_inca                                                           &
3467       &        (kdlon,kflev,dist, rmu0, fract, solaire,                             &
3468       &        paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,                 &
3469       &        wo(:, :, 1),                                                         &
3470       &        cldfrarad, cldemirad, cldtaurad,                                     &
3471       &        heat,heat0,cool,cool0,radsol,albpla,                                 &
3472       &        topsw,toplw,solsw,sollw,                                             &
3473       &        sollwdown,                                                           &
3474       &        topsw0,toplw0,solsw0,sollw0,                                         &
3475       &        lwdn0, lwdn, lwup0, lwup,                                            &
3476       &        swdn0, swdn, swup0, swup,                                            &
3477       &        ok_ade, ok_aie,                                                      &
3478       &        tau_aero, piz_aero, cg_aero,                                         &
3479       &        topswad_aero, solswad_aero,                                          &
3480       &        topswad0_aero, solswad0_aero,                                        &
3481       &        topsw_aero, topsw0_aero,                                             &
3482       &        solsw_aero, solsw0_aero,                                             &
3483       &        cldtaupirad,                                                         &
3484       &        topswai_aero, solswai_aero)
3485           
3486#endif
3487      ELSE
3488!c
3489!IM calcul radiatif pour le cas actuel
3490!c
3491       RCO2 = RCO2_act
3492       RCH4 = RCH4_act
3493       RN2O = RN2O_act
3494       RCFC11 = RCFC11_act
3495       RCFC12 = RCFC12_act
3496!c
3497      IF (prt_level .GE.10) THEN
3498       print *,' ->radlwsw, number 1 '
3499      ENDIF
3500!c
3501         CALL radlwsw                                                                &
3502       &        (dist, rmu0, fract,                                                  &
3503       &        paprs, pplay,zxtsol,albsol1, albsol2,                                &
3504       &        t_seri,q_seri,wo,                                                    &
3505       &        cldfrarad, cldemirad, cldtaurad,                                     &
3506       &        ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol,                  &
3507       &        flag_aerosol_strat,                                                  &
3508       &        tau_aero, piz_aero, cg_aero,                                         &
3509       &        cldtaupirad,new_aod,                                                 &
3510       &        zqsat, flwc, fiwc,                                                   &
3511       &        heat,heat0,cool,cool0,radsol,albpla,                                 &
3512       &        topsw,toplw,solsw,sollw,                                             &
3513       &        sollwdown,                                                           &
3514       &        topsw0,toplw0,solsw0,sollw0,                                         &
3515       &        lwdn0, lwdn, lwup0, lwup,                                            &
3516       &        swdn0, swdn, swup0, swup,                                            &
3517       &        topswad_aero, solswad_aero,                                          &
3518       &        topswai_aero, solswai_aero,                                          &
3519       &        topswad0_aero, solswad0_aero,                                        &
3520       &        topsw_aero, topsw0_aero,                                             &
3521       &        solsw_aero, solsw0_aero,                                             &
3522       &        topswcf_aero, solswcf_aero)
3523         
3524!c
3525!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3526!IM des taux doit etre different du taux actuel
3527!IM Par defaut on a les taux perturbes egaux aux taux actuels
3528!c
3529
3530      if (ok_4xCO2atm) then
3531       if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR.                          &
3532       &RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR.                         &
3533       &RCFC12_per.NE.RCFC12_act) THEN
3534!c
3535       RCO2 = RCO2_per
3536       RCH4 = RCH4_per
3537       RN2O = RN2O_per
3538       RCFC11 = RCFC11_per
3539       RCFC12 = RCFC12_per
3540!c
3541      IF (prt_level .GE.10) THEN
3542       print *,' ->radlwsw, number 2 '
3543      ENDIF
3544!c
3545         CALL radlwsw                                                                &
3546       &        (dist, rmu0, fract,                                                  &
3547       &        paprs, pplay,zxtsol,albsol1, albsol2,                                &
3548       &        t_seri,q_seri,wo,                                                    &
3549       &        cldfra, cldemi, cldtau,                                              &
3550       &        ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol,                  &
3551       &        flag_aerosol_strat,                                                  &
3552       &        tau_aero, piz_aero, cg_aero,                                         &
3553       &        cldtaupi,new_aod,                                                    &
3554       &        zqsat, flwc, fiwc,                                                   &
3555       &        heatp,heat0p,coolp,cool0p,radsolp,albplap,                           &
3556       &        topswp,toplwp,solswp,sollwp,                                         &
3557       &        sollwdownp,                                                          &
3558       &        topsw0p,toplw0p,solsw0p,sollw0p,                                     &
3559       &        lwdn0p, lwdnp, lwup0p, lwupp,                                        &
3560       &        swdn0p, swdnp, swup0p, swupp,                                        &
3561       &        topswad_aerop, solswad_aerop,                                        &
3562       &        topswai_aerop, solswai_aerop,                                        &
3563       &        topswad0_aerop, solswad0_aerop,                                      &
3564       &        topsw_aerop, topsw0_aerop,                                           &
3565       &        solsw_aerop, solsw0_aerop,                                           &
3566       &        topswcf_aerop, solswcf_aerop)
3567       endif
3568      endif
3569!c
3570      ENDIF ! aerosol_couple
3571      itaprad = 0
3572      ENDIF ! MOD(itaprad,radpas)
3573      itaprad = itaprad + 1
3574
3575      IF (iflag_radia.eq.0) THEN
3576         IF (prt_level.ge.9) THEN
3577            PRINT *,'--------------------------------------------------'
3578            PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3579            PRINT *,'>>>>           heat et cool mis a zero '
3580            PRINT *,'--------------------------------------------------'
3581         END IF
3582         heat=0.
3583         cool=0.
3584         sollw=0.   ! MPL 01032011
3585         solsw=0.
3586         radsol=0.
3587         swup=0.    ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
3588         swup0=0.
3589         swdn=0.
3590         swdn0=0.
3591         lwup=0.
3592         lwup0=0.
3593         lwdn=0.
3594         lwdn0=0.
3595      END IF
3596
3597!c
3598!c Ajouter la tendance des rayonnements (tous les pas)
3599!c
3600      DO k = 1, klev
3601      DO i = 1, klon
3602         t_seri(i,k) = t_seri(i,k)                                                   &
3603       &               + (heat(i,k)-cool(i,k)) * dtime/RDAY
3604      ENDDO
3605      ENDDO
3606!c
3607      if (mydebug) then
3608        call writefield_phy('u_seri',u_seri,llm)
3609        call writefield_phy('v_seri',v_seri,llm)
3610        call writefield_phy('t_seri',t_seri,llm)
3611        call writefield_phy('q_seri',q_seri,llm)
3612      endif
3613 
3614!IM
3615      IF (ip_ebil_phy.ge.2) THEN
3616        ztit='after rad'
3617        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime                             &
3618       &      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay              &
3619       &      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3620        call diagphy(airephy,ztit,ip_ebil_phy                                        &
3621       &      , topsw, toplw, solsw, sollw, zero_v                                   &
3622       &      , zero_v, zero_v, zero_v, ztsol                                        &
3623       &      , d_h_vcol, d_qt, d_ec                                                 &
3624       &      , fs_bound, fq_bound )
3625      END IF
3626      PRINT *,'  Lluis 9 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
3627!c
3628!c
3629!c Calculer l'hydrologie de la surface
3630!c
3631!c      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3632!c     .            agesno, ftsol,fqsurf,fsnow, ruis)
3633!c
3634
3635!c
3636!c Calculer le bilan du sol et la derive de temperature (couplage)
3637!c
3638      DO i = 1, klon
3639!c         bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3640!c a la demande de JLD
3641         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3642      ENDDO
3643!c
3644!cmoddeblott(jan95)
3645!c Appeler le programme de parametrisation de l'orographie
3646!c a l'echelle sous-maille:
3647!c
3648      IF (prt_level .GE.10) THEN
3649       print *,' call orography ? ', ok_orodr
3650      ENDIF
3651!c
3652      IF (ok_orodr) THEN
3653!c
3654!c  selection des points pour lesquels le shema est actif:
3655        igwd=0
3656        DO i=1,klon
3657        itest(i)=0
3658!c        IF ((zstd(i).gt.10.0)) THEN
3659        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3660          itest(i)=1
3661          igwd=igwd+1
3662          idx(igwd)=i
3663        ENDIF
3664        ENDDO
3665!c        igwdim=MAX(1,igwd)
3666!c
3667        IF (ok_strato) THEN
3668       
3669          CALL drag_noro_strato(klon,klev,dtime,paprs,pplay,                         &
3670       &                   zmea,zstd, zsig, zgam, zthe,zpic,zval,                    &
3671       &                   igwd,idx,itest,                                           &
3672       &                   t_seri, u_seri, v_seri,                                   &
3673       &                   zulow, zvlow, zustrdr, zvstrdr,                           &
3674       &                   d_t_oro, d_u_oro, d_v_oro)
3675
3676       ELSE
3677        CALL drag_noro(klon,klev,dtime,paprs,pplay,                                  &
3678       &                   zmea,zstd, zsig, zgam, zthe,zpic,zval,                    &
3679       &                   igwd,idx,itest,                                           &
3680       &                   t_seri, u_seri, v_seri,                                   &
3681       &                   zulow, zvlow, zustrdr, zvstrdr,                           &
3682       &                   d_t_oro, d_u_oro, d_v_oro)
3683       ENDIF
3684!c
3685!c  ajout des tendances
3686!-----------------------------------------------------------------------------------------
3687! ajout des tendances de la trainee de l'orographie
3688      CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
3689!-----------------------------------------------------------------------------------------
3690!c
3691      ENDIF ! fin de test sur ok_orodr
3692!c
3693      if (mydebug) then
3694        call writefield_phy('u_seri',u_seri,llm)
3695        call writefield_phy('v_seri',v_seri,llm)
3696        call writefield_phy('t_seri',t_seri,llm)
3697        call writefield_phy('q_seri',q_seri,llm)
3698      endif
3699     
3700      IF (ok_orolf) THEN
3701!c
3702!c  selection des points pour lesquels le shema est actif:
3703        igwd=0
3704        DO i=1,klon
3705        itest(i)=0
3706        IF ((zpic(i)-zmea(i)).GT.100.) THEN
3707          itest(i)=1
3708          igwd=igwd+1
3709          idx(igwd)=i
3710        ENDIF
3711        ENDDO
3712!c        igwdim=MAX(1,igwd)
3713!c
3714        IF (ok_strato) THEN
3715
3716          CALL lift_noro_strato(klon,klev,dtime,paprs,pplay,                         &
3717       &                   rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval,                  &
3718       &                   igwd,idx,itest,                                           &
3719       &                   t_seri, u_seri, v_seri,                                   &
3720       &                   zulow, zvlow, zustrli, zvstrli,                           &
3721       &                   d_t_lif, d_u_lif, d_v_lif               )
3722       
3723        ELSE
3724          CALL lift_noro(klon,klev,dtime,paprs,pplay,                                &
3725       &                   rlat,zmea,zstd,zpic,                                      &
3726       &                   itest,                                                    &
3727       &                   t_seri, u_seri, v_seri,                                   &
3728       &                   zulow, zvlow, zustrli, zvstrli,                           &
3729       &                   d_t_lif, d_u_lif, d_v_lif)
3730       ENDIF
3731!c   
3732!-----------------------------------------------------------------------------------------
3733! ajout des tendances de la portance de l'orographie
3734      CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3735!-----------------------------------------------------------------------------------------
3736!c
3737      ENDIF ! fin de test sur ok_orolf
3738!C  HINES GWD PARAMETRIZATION
3739
3740       IF (ok_hines) then
3741
3742         CALL hines_gwd(klon,klev,dtime,paprs,pplay,                                 &
3743       &                  rlat,t_seri,u_seri,v_seri,                                 &
3744       &                  zustrhi,zvstrhi,                                           &
3745       &                  d_t_hin, d_u_hin, d_v_hin)
3746!c
3747!c  ajout des tendances
3748        CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'hin')
3749
3750      ENDIF
3751!c
3752
3753!c
3754!IM cf. FLott BEG
3755!C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3756
3757      if (mydebug) then
3758        call writefield_phy('u_seri',u_seri,llm)
3759        call writefield_phy('v_seri',v_seri,llm)
3760        call writefield_phy('t_seri',t_seri,llm)
3761        call writefield_phy('q_seri',q_seri,llm)
3762      endif
3763
3764      DO i = 1, klon
3765        zustrph(i)=0.
3766        zvstrph(i)=0.
3767      ENDDO
3768      DO k = 1, klev
3769      DO i = 1, klon
3770       zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*                             &
3771       &            (paprs(i,k)-paprs(i,k+1))/rg
3772       zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*                             &
3773       &            (paprs(i,k)-paprs(i,k+1))/rg
3774      ENDDO
3775      ENDDO
3776!c
3777!IM calcul composantes axiales du moment angulaire et couple des montagnes
3778!c
3779!      IF (is_sequential .and. ok_orodr) THEN
3780!        CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,                            &
3781!       &                 ra,rg,romega,                                               &
3782!       &                 rlat,rlon,pphis,                                            &
3783!       &                 zustrdr,zustrli,zustrph,                                    &
3784!       &                 zvstrdr,zvstrli,zvstrph,                                    &
3785!       &                 paprs,u,v,                                                  &
3786!       &                 aam, torsfc)
3787!       ENDIF
3788!IM cf. FLott END
3789!IM
3790      IF (ip_ebil_phy.ge.2) THEN
3791        ztit='after orography'
3792        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime                             &
3793       &      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay              &
3794       &      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3795         call diagphy(airephy,ztit,ip_ebil_phy                                       &
3796       &      , zero_v, zero_v, zero_v, zero_v, zero_v                               &
3797       &      , zero_v, zero_v, zero_v, ztsol                                        &
3798       &      , d_h_vcol, d_qt, d_ec                                                 &
3799       &      , fs_bound, fq_bound )
3800      END IF
3801      PRINT *,'  Lluis 10 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
3802!c
3803!c
3804!====================================================================
3805! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
3806!====================================================================
3807! Abderrahmane 24.08.09
3808
3809      IF (ok_cosp) THEN
3810! adeclarer
3811#ifdef CPP_COSP
3812       IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
3813
3814       print*,'freq_cosp',freq_cosp
3815          mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
3816!       print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
3817!     s        ref_liq,ref_ice
3818          call phys_cosp(itap,dtime,freq_cosp,                                       &
3819       &                   ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP,                   &
3820       &                   ecrit_mth,ecrit_day,ecrit_hf,                             &
3821       &                   klon,klev,rlon,rlat,presnivs,overlap,                     &
3822       &                   ref_liq,ref_ice,                                          &
3823       &                   pctsrf(:,is_ter)+pctsrf(:,is_lic),                        &
3824       &                   zu10m,zv10m,pphis,                                        &
3825       &                   zphi,paprs(:,1:klev),pplay,zxtsol,t_seri,                 &
3826       &                   qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc,              &
3827       &                   prfl(:,1:klev),psfl(:,1:klev),                            &
3828       &                   pmflxr(:,1:klev),pmflxs(:,1:klev),                        &
3829       &                   mr_ozone,cldtau, cldemi)
3830
3831!     L          calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
3832!     L          cfaddbze,clcalipso2,dbze,cltlidarradar,
3833!     M          clMISR,
3834!     R          clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
3835!     I          tauisccp,albisccp,meantbisccp,meantbclrisccp)
3836
3837         ENDIF
3838
3839#endif
3840       ENDIF  !ok_cosp
3841!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3842!cAA
3843!cAA Installation de l'interface online-offline pour traceurs
3844!cAA
3845!c====================================================================
3846!c   Calcul  des tendances traceurs
3847!c====================================================================
3848!C
3849
3850       IF (type_trac=='repr') THEN
3851          sh_in(:,:) = q_seri(:,:)
3852       ELSE
3853          sh_in(:,:) = qx(:,:,ivap)
3854       END IF
3855
3856      call phytrac (                                                                 &
3857       &     itap,     days_elapsed+1,    jH_cur,   debut,                           &
3858       &     lafin,    dtime,     u, v,     t,                                       &
3859       &     paprs,    pplay,     pmfu,     pmfd,                                    &
3860       &     pen_u,    pde_u,     pen_d,    pde_d,                                   &
3861       &     cdragh,   coefh(:,:,is_ave),     fm_therm, entr_therm,                  &
3862       &     u1,       v1,        ftsol,    pctsrf,                                  &
3863       &     ustar,     u10m,      v10m,                                             &
3864       &     rlat,     rlon,                                                         &
3865       &     frac_impa,frac_nucl, beta_prec_fisrt,beta_prec,                         &
3866       &     presnivs, pphis,     pphi,     albsol1,                                 &
3867       &     sh_in,    rhcl,      cldfra,   rneb,                                    &
3868       &     diafra,   cldliq,    itop_con, ibas_con,                                &
3869       &     pmflxr,   pmflxs,    prfl,     psfl,                                    &
3870       &     da,       phi,       mp,       upwd,                                    &
3871       &     phi2,     d1a,       dam,      sij,                                     &
3872       &     wdtrainA, wdtrainM,  sigd,     clw,elij,                                &
3873       &     ev,       ep,        epmlmMm,  eplaMm,                                  &
3874       &     dnwd,     aerosol_couple,      flxmass_w,                               &
3875       &     tau_aero, piz_aero,  cg_aero,  ccm,                                     &
3876       &     rfname,                                                                 &
3877       &     d_tr_dyn,                                                               &
3878       &     tr_seri)
3879
3880      IF (offline) THEN
3881
3882       IF (prt_level.ge.9)                                                           &
3883       &    print*,'Attention on met a 0 les thermiques pour phystoke'
3884         call phystokenc (                                                                  &
3885       &                   nlon,klev,pdtphys,rlon,rlat,                              &
3886       &                   t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,                 &
3887       &                   fm_therm,entr_therm,                                      &
3888       &                   cdragh,coefh(:,:,is_ave),u1,v1,ftsol,pctsrf,              &
3889       &                   frac_impa, frac_nucl,                                     &
3890       &                   pphis,airephy,dtime,itap,                                 &
3891       &                   qx(:,:,ivap),da,phi,mp,upwd,dnwd)
3892
3893
3894      ENDIF
3895
3896!c
3897!c Calculer le transport de l'eau et de l'energie (diagnostique)
3898!c
3899      CALL transp (paprs,zxtsol,                                                     &
3900       &                   t_seri, q_seri, u_seri, v_seri, zphi,                     &
3901       &                   ve, vq, ue, uq)
3902!c
3903!IM global posePB BEG
3904      IF(1.EQ.0) THEN
3905!c
3906      CALL transp_lay (paprs,zxtsol,                                                 &
3907       &                   t_seri, q_seri, u_seri, v_seri, zphi,                     &
3908       &                   ve_lay, vq_lay, ue_lay, uq_lay)
3909!c
3910      ENDIF !(1.EQ.0) THEN
3911!IM global posePB END
3912!c Accumuler les variables a stocker dans les fichiers histoire:
3913!c
3914
3915!================================================================
3916! Conversion of kinetic and potential energy into heat, for
3917! parameterisation of subgrid-scale motions
3918!================================================================
3919
3920      d_t_ec(:,:)=0.
3921      forall (k=1: llm) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
3922      CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),                        &
3923       &        u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:),           &
3924       &        zmasse,exner,d_t_ec)
3925      t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
3926
3927!IM
3928      IF (ip_ebil_phy.ge.1) THEN
3929        ztit='after physic'
3930        CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime                             &
3931       &      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay              &
3932       &      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3933!C     Comme les tendances de la physique sont ajoute dans la dynamique,
3934!C     on devrait avoir que la variation d'entalpie par la dynamique
3935!C     est egale a la variation de la physique au pas de temps precedent.
3936!C     Donc la somme de ces 2 variations devrait etre nulle.
3937
3938       PRINT *,'  LLuis sending: d_h_vcol: ',d_h_vcol,' d_qt ',d_qt,' d_ec ',d_ec
3939
3940        call diagphy(airephy,ztit,ip_ebil_phy                                        &
3941       &      , topsw, toplw, solsw, sollw, sens                                     &
3942       &      , evap, rain_fall, snow_fall, ztsol                                    &
3943       &      , d_h_vcol, d_qt, d_ec                                                 &
3944       &      , fs_bound, fq_bound )
3945!C
3946      d_h_vcol_phy=d_h_vcol
3947!C
3948      END IF
3949      PRINT *,'  Lluis 11 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
3950
3951!C
3952!c=======================================================================
3953!c   SORTIES
3954!c=======================================================================
3955
3956!IM Interpolation sur les niveaux de pression du NMC
3957!c   -------------------------------------------------
3958!c
3959#include "calcul_STDlev.h"
3960!c
3961!c slp sea level pressure
3962      slp(:) = paprs(:,1)*exp(pphis(:)/(RD*t_seri(:,1)))
3963!c
3964!ccc prw = eau precipitable
3965      DO i = 1, klon
3966       prw(i) = 0.
3967       DO k = 1, klev
3968        prw(i) = prw(i) +                                                            &
3969       &           q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3970       ENDDO
3971      ENDDO
3972!c
3973!IM initialisation + calculs divers diag AMIP2
3974!c
3975#include "calcul_divers.h"
3976!c
3977      IF (type_trac == 'inca') THEN
3978#ifdef INCA
3979         CALL VTe(VTphysiq)
3980         CALL VTb(VTinca)
3981
3982         CALL chemhook_end (                                                         &
3983       &                        dtime,                                               &
3984       &                        pplay,                                               &
3985       &                        t_seri,                                              &
3986       &                        tr_seri,                                             &
3987       &                        nbtr,                                                &
3988       &                        paprs,                                               &
3989       &                        q_seri,                                              &
3990       &                        airephy,                                             &
3991       &                        pphi,                                                &
3992       &                        pphis,                                               &
3993       &                        zx_rh)
3994
3995         CALL VTe(VTinca)
3996         CALL VTb(VTphysiq)
3997#endif
3998      END IF
3999
4000
4001!c
4002!c Convertir les incrementations en tendances
4003!c
4004      IF (prt_level .GE.10) THEN
4005        print *,'Convertir les incrementations en tendances '
4006      ENDIF
4007!c
4008      if (mydebug) then
4009        call writefield_phy('u_seri',u_seri,llm)
4010        call writefield_phy('v_seri',v_seri,llm)
4011        call writefield_phy('t_seri',t_seri,llm)
4012        call writefield_phy('q_seri',q_seri,llm)
4013      endif
4014
4015      DO k = 1, klev
4016      DO i = 1, klon
4017         d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
4018         d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
4019         d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
4020         d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
4021         d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
4022      ENDDO
4023      ENDDO
4024!c
4025      IF (nqtot.GE.3) THEN
4026      DO iq = 3, nqtot
4027      DO  k = 1, klev
4028      DO  i = 1, klon
4029         d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
4030      ENDDO
4031      ENDDO
4032      ENDDO
4033      ENDIF
4034!c
4035!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
4036!IM global posePB#include "write_bilKP_ins.h"
4037!IM global posePB#include "write_bilKP_ave.h"
4038!c
4039
4040!c Sauvegarder les valeurs de t et q a la fin de la physique:
4041!c
4042      DO k = 1, klev
4043      DO i = 1, klon
4044         u_ancien(i,k) = u_seri(i,k)
4045         v_ancien(i,k) = v_seri(i,k)
4046         t_ancien(i,k) = t_seri(i,k)
4047         q_ancien(i,k) = q_seri(i,k)
4048      ENDDO
4049      ENDDO
4050
4051!!! RomP >>>
4052      IF (nqtot.GE.3) THEN
4053        DO iq = 3, nqtot
4054        DO k = 1, klev
4055        DO i = 1, klon
4056           tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
4057        ENDDO
4058        ENDDO
4059        ENDDO
4060      ENDIF
4061!!! RomP <<<
4062!==========================================================================
4063! Sorties des tendances pour un point particulier
4064! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
4065! pour le debug
4066! La valeur de igout est attribuee plus haut dans le programme
4067!==========================================================================
4068
4069      if (prt_level.ge.1) then
4070      write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
4071      write(lunout,*)                                                                &
4072       & 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
4073      write(lunout,*)                                                                &
4074       &  nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,                      &
4075       &  pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),           &
4076       &  pctsrf(igout,is_sic)
4077      write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
4078      do k=1,klev
4079         write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),                          &
4080       &   d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),                      &
4081       &   d_t_eva(igout,k)
4082      enddo
4083      write(lunout,*) 'cool,heat'
4084      do k=1,klev
4085         write(lunout,*) cool(igout,k),heat(igout,k)
4086      enddo
4087
4088      write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4089      do k=1,klev
4090         write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),                          &
4091       & d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4092      enddo
4093
4094      write(lunout,*) 'd_ps ',d_ps(igout)
4095      write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
4096      do k=1,klev,klon
4097         write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),                     &
4098       &  d_qx(igout,k,1),d_qx(igout,k,2)
4099      enddo
4100      endif
4101
4102!==========================================================================
4103
4104!c============================================================
4105!c   Calcul de la temperature potentielle
4106!c============================================================
4107      DO k = 1, klev
4108      DO i = 1, klon
4109!cJYG/IM theta en debut du pas de temps
4110!cJYG/IM       theta(i,k)=t(i,k)*(100000./pplay(i,k))**(RD/RCPD)
4111!cJYG/IM theta en fin de pas de temps de physique
4112        theta(i,k)=t_seri(i,k)*(100000./pplay(i,k))**(RD/RCPD)
4113!c thetal: 2 lignes suivantes a decommenter si vous avez les fichiers     MPL 20130625
4114!c fth_fonctions.F90 et parkind1.F90
4115!c sinon thetal=theta
4116!c       thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
4117!c    :         ql_seri(i,k))
4118        thetal(i,k)=theta(i,k)
4119      ENDDO
4120      ENDDO
4121!c
4122
4123!c 22.03.04 BEG
4124!c=============================================================
4125!c   Ecriture des sorties
4126!c=============================================================
4127#ifdef CPP_IOIPSL
4128 
4129!c Recupere des varibles calcule dans differents modules
4130!c pour ecriture dans histxxx.nc
4131
4132      ! Get some variables from module fonte_neige_mod
4133!L. Fita, LMD. November 2013
4134!! It is not working after removing starphy.nc and limit.nc?
4135      CALL fonte_neige_get_vars(pctsrf,                                              &
4136       &     zxfqcalving, zxfqfonte, zxffonte)
4137
4138
4139
4140!c=============================================================
4141! Separation entre thermiques et non thermiques dans les sorties
4142! de fisrtilp
4143!c=============================================================
4144
4145      if (iflag_thermals>=1) then
4146      d_t_lscth=0.
4147      d_t_lscst=0.
4148      d_q_lscth=0.
4149      d_q_lscst=0.
4150      do k=1,klev
4151         do i=1,klon
4152            if (ptconvth(i,k)) then
4153                d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4154                d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4155            else
4156                d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4157                d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4158            endif
4159         enddo
4160      enddo
4161
4162      do i=1,klon
4163      plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
4164      plul_th(i)=prfl(i,1)+psfl(i,1)
4165      enddo
4166      endif
4167
4168#include "phys_output_write_new.h"
4169
4170
4171#ifdef histISCCP
4172#include "write_histISCCP.h"
4173#endif
4174
4175#ifdef histNMC
4176#include "write_histhfNMC.h"
4177#include "write_histdayNMC.h"
4178#include "write_histmthNMC.h"
4179#endif
4180
4181#include "write_histday_seri.h"
4182
4183#include "write_paramLMDZ_phy.h"
4184
4185#endif
4186
4187!c 22.03.04 END
4188!c
4189!c====================================================================
4190!c Si c'est la fin, il faut conserver l'etat de redemarrage
4191!c====================================================================
4192!c
4193
4194!c        -----------------------------------------------------------------
4195!c        WSTATS: Saving statistics
4196!c        -----------------------------------------------------------------
4197!c        ("stats" stores and accumulates 8 key variables in file "stats.nc"
4198!c        which can later be used to make the statistic files of the run:
4199!c        "stats")          only possible in 3D runs !
4200
4201         
4202         IF (callstats) THEN
4203
4204           call wstats(klon,o_psol%name,"Surface pressure","Pa"                      &
4205       &                 ,2,paprs(:,1))
4206           call wstats(klon,o_tsol%name,"Surface temperature","K",                   &
4207       &                 2,zxtsol)
4208           zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
4209           call wstats(klon,o_precip%name,"Precip Totale liq+sol",                   &
4210       &                 "kg/(s*m2)",2,zx_tmp_fi2d)
4211           zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
4212           call wstats(klon,o_plul%name,"Large-scale Precip",                        &
4213       &                 "kg/(s*m2)",2,zx_tmp_fi2d)
4214           zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
4215           call wstats(klon,o_pluc%name,"Convective Precip",                         &
4216       &                 "kg/(s*m2)",2,zx_tmp_fi2d)
4217           call wstats(klon,o_sols%name,"Solar rad. at surf.",                       &
4218       &                 "W/m2",2,solsw)
4219           call wstats(klon,o_soll%name,"IR rad. at surf.",                          &
4220       &                 "W/m2",2,sollw)
4221          zx_tmp_fi2d(:) = topsw(:)-toplw(:)
4222          call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA",              &
4223       &                 "W/m2",2,zx_tmp_fi2d)
4224
4225
4226
4227           call wstats(klon,o_temp%name,"Air temperature","K",                       &
4228       &                 3,t_seri)
4229           call wstats(klon,o_vitu%name,"Zonal wind","m.s-1",                        &
4230       &                 3,u_seri)
4231           call wstats(klon,o_vitv%name,"Meridional wind",                           &
4232       &                "m.s-1",3,v_seri)
4233           call wstats(klon,o_vitw%name,"Vertical wind",                             &
4234       &                "m.s-1",3,omega)
4235           call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg",                &
4236       &                 3,q_seri)
4237 
4238
4239
4240           IF(lafin) THEN
4241             write (*,*) "Writing stats..."
4242             call mkstats(ierr)
4243           ENDIF
4244
4245         ENDIF !if callstats
4246
4247      IF (lafin) THEN
4248         itau_phy = itau_phy + itap
4249         CALL phyredem ("restartphy.nc")
4250!         open(97,form="unformatted",file="finbin")
4251!         write(97) u_seri,v_seri,t_seri,q_seri
4252!         close(97)
4253!$OMP MASTER
4254         if (read_climoz >= 1) then
4255            if (is_mpi_root) then
4256               call nf95_close(ncid_climoz)
4257            end if
4258            deallocate(press_climoz) ! pointer
4259         end if
4260!$OMP END MASTER
4261      ENDIF
4262     
4263!      first=.false.
4264
4265      RETURN
4266      END SUBROUTINE physiq
4267
4268      FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
4269      IMPLICIT none
4270!c
4271!c Calculer et imprimer l'eau totale. A utiliser pour verifier
4272!c la conservation de l'eau
4273!c
4274#include "YOMCST.h"
4275      INTEGER klon,klev
4276      REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
4277      REAL aire(klon)
4278      REAL qtotal, zx, qcheck
4279      INTEGER i, k
4280!c
4281      zx = 0.0
4282      DO i = 1, klon
4283         zx = zx + aire(i)
4284      ENDDO
4285      qtotal = 0.0
4286      DO k = 1, klev
4287      DO i = 1, klon
4288         qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)                                &
4289       &                     *(paprs(i,k)-paprs(i,k+1))/RG
4290      ENDDO
4291      ENDDO
4292!c
4293      qcheck = qtotal/zx
4294!c
4295      RETURN
4296      END
4297
4298      SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
4299      IMPLICIT none
4300!c
4301!c Tranformer une variable de la grille physique a
4302!c la grille d'ecriture
4303!c
4304      INTEGER nfield,nlon,iim,jjmp1, jjm
4305      REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
4306!c
4307      INTEGER i, n, ig
4308!c
4309      jjm = jjmp1 - 1
4310      DO n = 1, nfield
4311         DO i=1,iim
4312            ecrit(i,n) = fi(1,n)
4313            ecrit(i+jjm*iim,n) = fi(nlon,n)
4314         ENDDO
4315         DO ig = 1, nlon - 2
4316           ecrit(iim+ig,n) = fi(1+ig,n)
4317         ENDDO
4318      ENDDO
4319      RETURN
4320      END SUBROUTINE gr_fi_ecrit
Note: See TracBrowser for help on using the repository browser.