source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/physiq_F90_v2039 @ 2404

Last change on this file since 2404 was 2175, checked in by jescribano, 10 years ago

SPLA code included for first time

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