source: LMDZ5/trunk/libf/phylmd/physiq.F90 @ 2003

Last change on this file since 2003 was 2003, checked in by Laurent Fairhead, 10 years ago

Nouvelle version qui inclut les effets des aérosols et propose les mêmes diagnostics des effets
directs et indirects que l'ancienne version du rayonnement.
OB

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