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

Last change on this file since 2006 was 2006, checked in by fhourdin, 10 years ago

Modification de la spécification de la plage de température pour
la phase mixte liquide/glace des nuages.
Contrôle par les paramètres t_glace_min/max, exposant_glace, iflag_t_glace

Modifying the specification of the mixte liquid/ice phase for cloud water.

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