source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/physiq.F90 @ 2404

Last change on this file since 2404 was 2303, checked in by jescribano, 9 years ago

Bugs corrections, control vector is now fine mode+coarse mode and seasalt coarse+fine, change in emission scheme parameters, more outputs at 10h30 and 13h30 LT. (Pending correct optical and sedimentation parameters)

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