source: LMDZ5/trunk/libf/phylmd/physiq.F @ 1454

Last change on this file since 1454 was 1454, checked in by Laurent Fairhead, 14 years ago

Merge of LMDZ5V1.0-dev branch r1453 into LMDZ5 trunk r1434


Fusion entre la version r1453 de la branche de développement LMDZ5V1.0-dev
et le tronc LMDZ5 (r1434)

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