source: LMDZ4/trunk/libf/phylmd/physiq.F @ 1645

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

Problems with some initialisations in the new physics package. These are
temporary fixes. Further investigations needed. Trac ticket issued. JYG


Problèmes d'initialisations dans la nouvelle physique. Ces corrections
sont temporaires, il faut trouver la vraie raison des plantages. Un ticket
trac est créé. JYG

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 123.0 KB
Line 
1! $Id: physiq.F 1428 2010-09-13 08:43:37Z lguez $
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
1252          paire_ter(:)=0.   
1253cIM 180608
1254c         pmflxr=0.
1255c         pmflxs=0.
1256
1257        itau_con=0
1258        first=.false.
1259
1260      endif  ! first
1261
1262       modname = 'physiq'
1263cIM
1264      IF (ip_ebil_phy.ge.1) THEN
1265        DO i=1,klon
1266          zero_v(i)=0.
1267        END DO
1268      END IF
1269      ok_sync=.TRUE.
1270
1271      IF (debut) THEN
1272         CALL suphel ! initialiser constantes et parametres phys.
1273      ENDIF
1274
1275      if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
1276
1277
1278c======================================================================
1279! Gestion calendrier : mise a jour du module phys_cal_mod
1280!
1281c     CALL phys_cal_update(jD_cur,jH_cur)
1282
1283c
1284c Si c'est le debut, il faut initialiser plusieurs choses
1285c          ********
1286c
1287       IF (debut) THEN
1288!rv
1289cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
1290cde la convection a partir des caracteristiques du thermique
1291         wght_th(:,:)=1.
1292         lalim_conv(:)=1
1293cRC
1294         u10m(:,:)=0.
1295         v10m(:,:)=0.
1296         rain_con(:)=0.
1297         snow_con(:)=0.
1298         topswai(:)=0.
1299         topswad(:)=0.
1300         solswai(:)=0.
1301         solswad(:)=0.
1302
1303         lambda_th(:,:)=0.
1304         wmax_th(:)=0.
1305         tau_overturning_th(:)=0.
1306
1307         IF (config_inca /= 'none') THEN
1308            ! jg : initialisation jusqu'au ces variables sont dans restart
1309            ccm(:,:,:) = 0.
1310            tau_aero(:,:,:,:) = 0.
1311            piz_aero(:,:,:,:) = 0.
1312            cg_aero(:,:,:,:) = 0.
1313         END IF
1314
1315         rnebcon0(:,:) = 0.0
1316         clwcon0(:,:) = 0.0
1317         rnebcon(:,:) = 0.0
1318         clwcon(:,:) = 0.0
1319
1320cIM     
1321         IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1322c
1323      print*,'iflag_coupl,iflag_clos,iflag_wake',
1324     .   iflag_coupl,iflag_clos,iflag_wake
1325      print*,'CYCLE_DIURNE', cycle_diurne
1326c
1327      IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
1328         abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
1329         CALL abort_gcm (modname,abort_message,1)
1330      ENDIF
1331c
1332      IF(ok_isccp.AND.iflag_con.LE.2) THEN
1333         abort_message = 'ISCCP-like outputs may be available for KE
1334     .(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
1335         CALL abort_gcm (modname,abort_message,1)
1336      ENDIF
1337c
1338c Initialiser les compteurs:
1339c
1340         itap    = 0
1341         itaprad = 0
1342
1343!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1344!! Un petit travail à faire ici.
1345!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1346
1347         if (iflag_pbl>1) then
1348             PRINT*, "Using method MELLOR&YAMADA"
1349         endif
1350
1351!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1352! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1353! dyn3d
1354! Attention : la version precedente n'etait pas tres propre.
1355! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1356! pour obtenir le meme resultat.
1357         dtime=pdtphys
1358         radpas = NINT( 86400./dtime/nbapp_rad)
1359!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1360
1361         CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1362cIM begin
1363          print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)
1364     $,ratqs(1,1)
1365cIM end
1366
1367
1368
1369!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1370c
1371C on remet le calendrier a zero
1372c
1373         IF (raz_date .eq. 1) THEN
1374           itau_phy = 0
1375         ENDIF
1376
1377cIM cf. AM 081204 BEG
1378         PRINT*,'cycle_diurne3 =',cycle_diurne
1379cIM cf. AM 081204 END
1380c
1381         CALL printflag( tabcntr0,radpas,ok_journe,
1382     ,                    ok_instan, ok_region )
1383c
1384         IF (ABS(dtime-pdtphys).GT.0.001) THEN
1385            WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1386     .                        pdtphys
1387            abort_message='Pas physique n est pas correct '
1388!           call abort_gcm(modname,abort_message,1)
1389            dtime=pdtphys
1390         ENDIF
1391         IF (nlon .NE. klon) THEN
1392            WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1393     .                      klon
1394            abort_message='nlon et klon ne sont pas coherents'
1395            call abort_gcm(modname,abort_message,1)
1396         ENDIF
1397         IF (nlev .NE. klev) THEN
1398            WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1399     .                       klev
1400            abort_message='nlev et klev ne sont pas coherents'
1401            call abort_gcm(modname,abort_message,1)
1402         ENDIF
1403c
1404         IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
1405           WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1406           WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1407           abort_message='Nbre d appels au rayonnement insuffisant'
1408           call abort_gcm(modname,abort_message,1)
1409         ENDIF
1410         WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1411         WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1412     .                   ok_cvl
1413c
1414cKE43
1415c Initialisation pour la convection de K.E. (sb):
1416         IF (iflag_con.GE.3) THEN
1417
1418         WRITE(lunout,*)"*** Convection de Kerry Emanuel 4.3  "
1419         WRITE(lunout,*)
1420     .      "On va utiliser le melange convectif des traceurs qui"
1421         WRITE(lunout,*)"est calcule dans convect4.3"
1422         WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1423
1424          DO i = 1, klon
1425           ema_cbmf(i) = 0.
1426           ema_pcb(i)  = 0.
1427           ema_pct(i)  = 0.
1428c          ema_workcbmf(i) = 0.
1429          ENDDO
1430cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1431          DO i = 1, klon
1432           ibas_con(i) = 1
1433           itop_con(i) = 1
1434          ENDDO
1435cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1436c===============================================================================
1437cCR:04.12.07: initialisations poches froides
1438c Controle de ALE et ALP pour la fermeture convective (jyg)
1439          if (iflag_wake.eq.1) then
1440            CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
1441     s                  ,alp_bl_prescr, ale_bl_prescr)
1442c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1443c        print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1444          endif
1445
1446        do i = 1,klon
1447         Ale_bl(i)=0.
1448         Alp_bl(i)=0.
1449        enddo
1450
1451c================================================================================
1452
1453         ENDIF !debut
1454
1455           DO i=1,klon
1456             rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1457           ENDDO
1458
1459c34EK
1460         IF (ok_orodr) THEN
1461
1462!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1463! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1464! justement quand ok_orodr = false.
1465! ce rugoro est utilise par la couche limite et fait double emploi
1466! avec les paramétrisations spécifiques de Francois Lott.
1467!           DO i=1,klon
1468!             rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1469!           ENDDO
1470!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1471           IF (ok_strato) THEN
1472             CALL SUGWD_strato(klon,klev,paprs,pplay)
1473           ELSE
1474             CALL SUGWD(klon,klev,paprs,pplay)
1475           ENDIF
1476           
1477           DO i=1,klon
1478             zuthe(i)=0.
1479             zvthe(i)=0.
1480             if(zstd(i).gt.10.)then
1481               zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1482               zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1483             endif
1484           ENDDO
1485         ENDIF
1486c
1487c
1488         lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
1489         WRITE(lunout,*)'La frequence de lecture surface est de ',
1490     .                   lmt_pas
1491c
1492cIM 030306 END
1493
1494      capemaxcels = 't_max(X)'
1495      t2mincels = 't_min(X)'
1496      t2maxcels = 't_max(X)'
1497      tinst = 'inst(X)'
1498      tave = 'ave(X)'
1499cIM cf. AM 081204 BEG
1500      write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1501cIM cf. AM 081204 END
1502c
1503c=============================================================
1504c   Initialisation des sorties
1505c=============================================================
1506
1507#ifdef CPP_IOIPSL
1508
1509c$OMP MASTER
1510       call phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta,
1511     &                        ctetaSTD,dtime,ok_veget,
1512     &                        type_ocean,iflag_pbl,ok_mensuel,ok_journe,
1513     &                        ok_hf,ok_instan,ok_LES,ok_ade,ok_aie,
1514     &                        read_climoz, new_aod, aerosol_couple
1515     &                        )
1516c$OMP END MASTER
1517c$OMP BARRIER
1518
1519#ifdef histISCCP
1520#include "ini_histISCCP.h"
1521#endif
1522
1523#ifdef histNMC
1524#include "ini_histhfNMC.h"
1525#include "ini_histdayNMC.h"
1526#include "ini_histmthNMC.h"
1527#endif
1528
1529#include "ini_histday_seri.h"
1530
1531#include "ini_paramLMDZ_phy.h"
1532
1533#endif
1534
1535cIM 250308bad guide        ecrit_hf2mth = 30*1/ecrit_hf
1536         ecrit_hf2mth = ecrit_mth/ecrit_hf
1537
1538         ecrit_hf = ecrit_hf * un_jour
1539cIM
1540         IF(ecrit_day.LE.1.) THEN
1541          ecrit_day = ecrit_day * un_jour !en secondes
1542         ENDIF
1543cIM
1544         ecrit_mth = ecrit_mth * un_jour
1545         ecrit_ins = ecrit_ins * un_jour
1546         ecrit_reg = ecrit_reg * un_jour
1547         ecrit_tra = ecrit_tra * un_jour
1548         ecrit_ISCCP = ecrit_ISCCP * un_jour
1549         ecrit_LES = ecrit_LES * un_jour
1550c
1551         PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
1552     .   ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
1553     .   ecrit_hf2mth
1554cIM 030306 END
1555
1556
1557cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1558      date0 = jD_ref
1559      WRITE(*,*) 'physiq date0 : ',date0
1560c
1561c
1562c
1563c Prescrire l'ozone dans l'atmosphere
1564c
1565c
1566cc         DO i = 1, klon
1567cc         DO k = 1, klev
1568cc            CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1569cc         ENDDO
1570cc         ENDDO
1571c
1572      IF (config_inca /= 'none') THEN
1573#ifdef INCA
1574         CALL VTe(VTphysiq)
1575         CALL VTb(VTinca)
1576!         iii = MOD(NINT(xjour),360)
1577!         calday = REAL(iii) + jH_cur
1578         calday = REAL(days_elapsed) + jH_cur
1579         WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1580
1581         CALL chemini(
1582     $                   rg,
1583     $                   ra,
1584     $                   airephy,
1585     $                   rlat,
1586     $                   rlon,
1587     $                   presnivs,
1588     $                   calday,
1589     $                   klon,
1590     $                   nqtot,
1591     $                   pdtphys,
1592     $                   annee_ref,
1593     $                   day_ref,
1594     $                   itau_phy)
1595
1596         CALL VTe(VTinca)
1597         CALL VTb(VTphysiq)
1598#endif
1599      END IF
1600c
1601c
1602!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1603! Nouvelle initialisation pour le rayonnement RRTM
1604!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1605
1606      call iniradia(klon,klev,paprs(1,1:klev+1))
1607
1608C$omp single
1609      if (read_climoz >= 1) then
1610         call open_climoz(ncid_climoz, press_climoz)
1611      END IF
1612C$omp end single
1613      ENDIF
1614!
1615!   ****************     Fin  de   IF ( debut  )   ***************
1616!
1617!
1618! Incrementer le compteur de la physique
1619!
1620      itap   = itap + 1
1621!
1622! Update fraction of the sub-surfaces (pctsrf) and
1623! initialize, where a new fraction has appeared, all variables depending
1624! on the surface fraction.
1625!
1626      CALL change_srf_frac(itap, dtime, days_elapsed+1,
1627     *     pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
1628
1629! Tendances bidons pour les processus qui n'affectent pas certaines
1630! variables.
1631      du0(:,:)=0.
1632      dv0(:,:)=0.
1633      dq0(:,:)=0.
1634      dql0(:,:)=0.
1635c
1636c Mettre a zero des variables de sortie (pour securite)
1637c
1638      DO i = 1, klon
1639         d_ps(i) = 0.0
1640      ENDDO
1641      DO k = 1, klev
1642      DO i = 1, klon
1643         d_t(i,k) = 0.0
1644         d_u(i,k) = 0.0
1645         d_v(i,k) = 0.0
1646      ENDDO
1647      ENDDO
1648      DO iq = 1, nqtot
1649      DO k = 1, klev
1650      DO i = 1, klon
1651         d_qx(i,k,iq) = 0.0
1652      ENDDO
1653      ENDDO
1654      ENDDO
1655      da(:,:)=0.
1656      mp(:,:)=0.
1657      phi(:,:,:)=0.
1658c
1659c Ne pas affecter les valeurs entrees de u, v, h, et q
1660c
1661      DO k = 1, klev
1662      DO i = 1, klon
1663         t_seri(i,k)  = t(i,k)
1664         u_seri(i,k)  = u(i,k)
1665         v_seri(i,k)  = v(i,k)
1666         q_seri(i,k)  = qx(i,k,ivap)
1667         ql_seri(i,k) = qx(i,k,iliq)
1668         qs_seri(i,k) = 0.
1669      ENDDO
1670      ENDDO
1671      IF (nqtot.GE.3) THEN
1672      DO iq = 3, nqtot
1673      DO  k = 1, klev
1674      DO  i = 1, klon
1675         tr_seri(i,k,iq-2) = qx(i,k,iq)
1676      ENDDO
1677      ENDDO
1678      ENDDO
1679      ELSE
1680      DO k = 1, klev
1681      DO i = 1, klon
1682         tr_seri(i,k,1) = 0.0
1683      ENDDO
1684      ENDDO
1685      ENDIF
1686C
1687      DO i = 1, klon
1688        ztsol(i) = 0.
1689      ENDDO
1690      DO nsrf = 1, nbsrf
1691        DO i = 1, klon
1692          ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1693        ENDDO
1694      ENDDO
1695cIM
1696      IF (ip_ebil_phy.ge.1) THEN
1697        ztit='after dynamic'
1698        CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1699     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1700     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1701C     Comme les tendances de la physique sont ajoute dans la dynamique,
1702C     on devrait avoir que la variation d'entalpie par la dynamique
1703C     est egale a la variation de la physique au pas de temps precedent.
1704C     Donc la somme de ces 2 variations devrait etre nulle.
1705        call diagphy(airephy,ztit,ip_ebil_phy
1706     e      , zero_v, zero_v, zero_v, zero_v, zero_v
1707     e      , zero_v, zero_v, zero_v, ztsol
1708     e      , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1709     s      , fs_bound, fq_bound )
1710      END IF
1711
1712c Diagnostiquer la tendance dynamique
1713c
1714      IF (ancien_ok) THEN
1715         DO k = 1, klev
1716         DO i = 1, klon
1717            d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1718            d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1719            d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1720            d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1721         ENDDO
1722         ENDDO
1723      ELSE
1724         DO k = 1, klev
1725         DO i = 1, klon
1726            d_u_dyn(i,k) = 0.0
1727            d_v_dyn(i,k) = 0.0
1728            d_t_dyn(i,k) = 0.0
1729            d_q_dyn(i,k) = 0.0
1730         ENDDO
1731         ENDDO
1732         ancien_ok = .TRUE.
1733      ENDIF
1734c
1735c Ajouter le geopotentiel du sol:
1736c
1737      DO k = 1, klev
1738      DO i = 1, klon
1739         zphi(i,k) = pphi(i,k) + pphis(i)
1740      ENDDO
1741      ENDDO
1742c
1743c Verifier les temperatures
1744c
1745cIM BEG
1746      IF (check) THEN
1747       amn=MIN(ftsol(1,is_ter),1000.)
1748       amx=MAX(ftsol(1,is_ter),-1000.)
1749       DO i=2, klon
1750        amn=MIN(ftsol(i,is_ter),amn)
1751        amx=MAX(ftsol(i,is_ter),amx)
1752       ENDDO
1753c
1754       PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1755      ENDIF !(check) THEN
1756cIM END
1757c
1758      CALL hgardfou(t_seri,ftsol,'debutphy')
1759c
1760cIM BEG
1761      IF (check) THEN
1762       amn=MIN(ftsol(1,is_ter),1000.)
1763       amx=MAX(ftsol(1,is_ter),-1000.)
1764       DO i=2, klon
1765        amn=MIN(ftsol(i,is_ter),amn)
1766        amx=MAX(ftsol(i,is_ter),amx)
1767       ENDDO
1768c
1769       PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1770      ENDIF !(check) THEN
1771cIM END
1772c
1773c Mettre en action les conditions aux limites (albedo, sst, etc.).
1774c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1775c
1776      if (read_climoz >= 1) then
1777C        Ozone from a file
1778!        Update required ozone index:
1779         ro3i = int((days_elapsed + jh_cur - jh_1jan)
1780     $        / ioget_year_len(year_cur) * 360.) + 1
1781         if (ro3i == 361) ro3i = 360
1782C        (This should never occur, except perhaps because of roundup
1783C        error. See documentation.)
1784         if (ro3i /= co3i) then
1785C           Update ozone field:
1786            if (read_climoz == 1) then
1787               call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i,
1788     $              press_in_edg=press_climoz, paprs=paprs, v3=wo)
1789            else
1790C              read_climoz == 2
1791               call regr_pr_av(ncid_climoz,
1792     $              (/"tro3         ", "tro3_daylight"/),
1793     $              julien=ro3i, press_in_edg=press_climoz, paprs=paprs,
1794     $              v3=wo)
1795            end if
1796!           Convert from mole fraction of ozone to column density of ozone in a
1797!           cell, in kDU:
1798            forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l)
1799     $           * rmo3 / rmd * zmasse / dobson_u / 1e3
1800C           (By regridding ozone values for LMDZ only once every 360th of
1801C           year, we have already neglected the variation of pressure in one
1802C           360th of year. So do not recompute "wo" at each time step even if
1803C           "zmasse" changes a little.)
1804            co3i = ro3i
1805         end if
1806      elseif (MOD(itap-1,lmt_pas) == 0) THEN
1807C        Once per day, update ozone from Royer:
1808         wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
1809      ENDIF
1810c
1811c Re-evaporer l'eau liquide nuageuse
1812c
1813      DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
1814      DO i = 1, klon
1815         zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1816c        zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1817         zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1818         zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1819         zb = MAX(0.0,ql_seri(i,k))
1820         za = - MAX(0.0,ql_seri(i,k))
1821     .                  * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
1822         t_seri(i,k) = t_seri(i,k) + za
1823         q_seri(i,k) = q_seri(i,k) + zb
1824         ql_seri(i,k) = 0.0
1825         d_t_eva(i,k) = za
1826         d_q_eva(i,k) = zb
1827      ENDDO
1828      ENDDO
1829cIM
1830      IF (ip_ebil_phy.ge.2) THEN
1831        ztit='after reevap'
1832        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1833     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1834     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1835         call diagphy(airephy,ztit,ip_ebil_phy
1836     e      , zero_v, zero_v, zero_v, zero_v, zero_v
1837     e      , zero_v, zero_v, zero_v, ztsol
1838     e      , d_h_vcol, d_qt, d_ec
1839     s      , fs_bound, fq_bound )
1840C
1841      END IF
1842
1843c
1844c=========================================================================
1845! Calculs de l'orbite.
1846! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1847! doit donc etre placé avant radlwsw et pbl_surface
1848
1849! calcul selon la routine utilisee pour les planetes
1850      if (new_orbit) then
1851        call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
1852        day_since_equinox = (jD_cur + jH_cur) - jD_eq
1853!        day_since_equinox = (jD_cur) - jD_eq
1854        call solarlong(day_since_equinox, zlongi, dist)
1855      else     
1856! calcul selon la routine utilisee pour l'AR4
1857!   choix entre calcul de la longitude solaire vraie ou valeur fixee a
1858!   solarlong0
1859        if (solarlong0<-999.) then
1860           CALL orbite(REAL(days_elapsed+1),zlongi,dist)
1861        else
1862           zlongi=solarlong0  ! longitude solaire vraie
1863           dist=1.            ! distance au soleil / moyenne
1864        endif
1865      endif
1866      if(prt_level.ge.1)                                                &
1867     &    write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
1868
1869!  Avec ou sans cycle diurne
1870      IF (cycle_diurne) THEN
1871        zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1872        CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
1873      ELSE
1874        CALL angle(zlongi, rlat, fract, rmu0)
1875      ENDIF
1876
1877      if (mydebug) then
1878        call writefield_phy('u_seri',u_seri,llm)
1879        call writefield_phy('v_seri',v_seri,llm)
1880        call writefield_phy('t_seri',t_seri,llm)
1881        call writefield_phy('q_seri',q_seri,llm)
1882      endif
1883
1884ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1885c Appel au pbl_surface : Planetary Boudary Layer et Surface
1886c Cela implique tous les interactions des sous-surfaces et la partie diffusion
1887c turbulent du couche limit.
1888c
1889c Certains varibales de sorties de pbl_surface sont utiliser que pour
1890c ecriture des fihiers hist_XXXX.nc, ces sont :
1891c   qsol,      zq2m,      s_pblh,  s_lcl,
1892c   s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
1893c   s_therm,   s_trmb1,   s_trmb2, s_trmb3,
1894c   zxrugs,    zu10m,     zv10m,   fder,
1895c   zxqsurf,   rh2m,      zxfluxu, zxfluxv,
1896c   frugs,     agesno,    fsollw,  fsolsw,
1897c   d_ts,      fevap,     fluxlat, t2m,
1898c   wfbils,    wfbilo,    fluxt,   fluxu, fluxv,
1899c
1900c Certains ne sont pas utiliser du tout :
1901c   dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1902c
1903
1904      CALL pbl_surface(
1905     e     dtime,     date0,     itap,    days_elapsed+1,
1906     e     debut,     lafin,
1907     e     rlon,      rlat,      rugoro,  rmu0,     
1908     e     rain_fall, snow_fall, solsw,   sollw,   
1909     e     t_seri,    q_seri,    u_seri,  v_seri,   
1910     e     pplay,     paprs,     pctsrf,           
1911     +     ftsol,     falb1,     falb2,   u10m,   v10m,
1912     s     sollwdown, cdragh,    cdragm,  u1,    v1,
1913     s     albsol1,   albsol2,   sens,    evap, 
1914     s     zxtsol,    zxfluxlat, zt2m,    qsat2m,
1915     s     d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf,
1916     s     coefh,     slab_wfbils,               
1917     d     qsol,      zq2m,      s_pblh,  s_lcl,
1918     d     s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
1919     d     s_therm,   s_trmb1,   s_trmb2, s_trmb3,
1920     d     zxrugs,    zu10m,     zv10m,   fder,
1921     d     zxqsurf,   rh2m,      zxfluxu, zxfluxv,
1922     d     frugs,     agesno,    fsollw,  fsolsw,
1923     d     d_ts,      fevap,     fluxlat, t2m,
1924     d     wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,
1925     -     dsens,     devap,     zxsnow,
1926     -     zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
1927
1928
1929!-----------------------------------------------------------------------------------------
1930! ajout des tendances de la diffusion turbulente
1931      CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
1932!-----------------------------------------------------------------------------------------
1933
1934      if (mydebug) then
1935        call writefield_phy('u_seri',u_seri,llm)
1936        call writefield_phy('v_seri',v_seri,llm)
1937        call writefield_phy('t_seri',t_seri,llm)
1938        call writefield_phy('q_seri',q_seri,llm)
1939      endif
1940
1941
1942      IF (ip_ebil_phy.ge.2) THEN
1943        ztit='after surface_main'
1944        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
1945     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1946     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1947         call diagphy(airephy,ztit,ip_ebil_phy
1948     e      , zero_v, zero_v, zero_v, zero_v, sens
1949     e      , evap  , zero_v, zero_v, ztsol
1950     e      , d_h_vcol, d_qt, d_ec
1951     s      , fs_bound, fq_bound )
1952      END IF
1953
1954c =================================================================== c
1955c   Calcul de Qsat
1956
1957      DO k = 1, klev
1958      DO i = 1, klon
1959         zx_t = t_seri(i,k)
1960         IF (thermcep) THEN
1961            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1962            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1963            zx_qs  = MIN(0.5,zx_qs)
1964            zcor   = 1./(1.-retv*zx_qs)
1965            zx_qs  = zx_qs*zcor
1966         ELSE
1967           IF (zx_t.LT.t_coup) THEN
1968              zx_qs = qsats(zx_t)/pplay(i,k)
1969           ELSE
1970              zx_qs = qsatl(zx_t)/pplay(i,k)
1971           ENDIF
1972         ENDIF
1973         zqsat(i,k)=zx_qs
1974      ENDDO
1975      ENDDO
1976
1977      if (prt_level.ge.1) then
1978      write(lunout,*) 'L   qsat (g/kg) avant clouds_gno'
1979      write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1980      endif
1981c
1982c Appeler la convection (au choix)
1983c
1984      DO k = 1, klev
1985      DO i = 1, klon
1986         conv_q(i,k) = d_q_dyn(i,k)
1987     .               + d_q_vdf(i,k)/dtime
1988         conv_t(i,k) = d_t_dyn(i,k)
1989     .               + d_t_vdf(i,k)/dtime
1990      ENDDO
1991      ENDDO
1992      IF (check) THEN
1993         za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1994         WRITE(lunout,*) "avantcon=", za
1995      ENDIF
1996      zx_ajustq = .FALSE.
1997      IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1998      IF (zx_ajustq) THEN
1999         DO i = 1, klon
2000            z_avant(i) = 0.0
2001         ENDDO
2002         DO k = 1, klev
2003         DO i = 1, klon
2004            z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
2005     .                        *(paprs(i,k)-paprs(i,k+1))/RG
2006         ENDDO
2007         ENDDO
2008      ENDIF
2009
2010c Calcule de vitesse verticale a partir de flux de masse verticale
2011      DO k = 1, klev
2012         DO i = 1, klon
2013            omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
2014         END DO
2015      END DO
2016      if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ',
2017     $     omega(igout, :)
2018
2019      IF (iflag_con.EQ.1) THEN
2020        abort_message ='reactiver le call conlmd dans physiq.F'
2021        CALL abort_gcm (modname,abort_message,1)
2022c     CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2023c    .             d_t_con, d_q_con,
2024c    .             rain_con, snow_con, ibas_con, itop_con)
2025      ELSE IF (iflag_con.EQ.2) THEN
2026      CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
2027     e            conv_t, conv_q, -evap, omega,
2028     s            d_t_con, d_q_con, rain_con, snow_con,
2029     s            pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2030     s            kcbot, kctop, kdtop, pmflxr, pmflxs)
2031      d_u_con = 0.
2032      d_v_con = 0.
2033
2034      WHERE (rain_con < 0.) rain_con = 0.
2035      WHERE (snow_con < 0.) snow_con = 0.
2036      DO i = 1, klon
2037         ibas_con(i) = klev+1 - kcbot(i)
2038         itop_con(i) = klev+1 - kctop(i)
2039      ENDDO
2040      ELSE IF (iflag_con.GE.3) THEN
2041c nb of tracers for the KE convection:
2042c MAF la partie traceurs est faite dans phytrac
2043c on met ntra=1 pour limiter les appels mais on peut
2044c supprimer les calculs / ftra.
2045              ntra = 1
2046
2047c=====================================================================================
2048cajout pour la parametrisation des poches froides:
2049ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
2050      do k=1,klev
2051            do i=1,klon
2052             if (iflag_wake.eq.1) then
2053             t_wake(i,k) = t_seri(i,k)
2054     .           +(1-wake_s(i))*wake_deltat(i,k)
2055             q_wake(i,k) = q_seri(i,k)
2056     .           +(1-wake_s(i))*wake_deltaq(i,k)
2057             t_undi(i,k) = t_seri(i,k)
2058     .           -wake_s(i)*wake_deltat(i,k)
2059             q_undi(i,k) = q_seri(i,k)
2060     .           -wake_s(i)*wake_deltaq(i,k)
2061             else
2062             t_wake(i,k) = t_seri(i,k)
2063             q_wake(i,k) = q_seri(i,k)
2064             t_undi(i,k) = t_seri(i,k)
2065             q_undi(i,k) = q_seri(i,k)
2066             endif
2067            enddo
2068         enddo
2069     
2070cc--   Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
2071cc--    pour le soulevement des particules dans le modele convectif
2072c
2073      do i = 1,klon
2074        ALE(i) = 0.
2075        ALP(i) = 0.
2076      enddo
2077c
2078ccalcul de ale_wake et alp_wake
2079       if (iflag_wake.eq.1) then
2080         if (itap .le. it_wape_prescr) then
2081          do i = 1,klon
2082           ale_wake(i) = wape_prescr
2083           alp_wake(i) = fip_prescr
2084          enddo
2085         else
2086          do i = 1,klon
2087cjyg  ALE=WAPE au lieu de ALE = 1/2 Cstar**2
2088ccc           ale_wake(i) = 0.5*wake_cstar(i)**2
2089           ale_wake(i) = wake_pe(i)
2090           alp_wake(i) = wake_fip(i)
2091          enddo
2092         endif
2093       else
2094         do i = 1,klon
2095           ale_wake(i) = 0.
2096           alp_wake(i) = 0.
2097         enddo
2098       endif
2099ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2100cdans le thermique sinon
2101       if (iflag_coupl.eq.0) then
2102          if (debut.and.prt_level.gt.9)
2103     $                     WRITE(lunout,*)'ALE et ALP imposes'
2104          do i = 1,klon
2105con ne couple que ale
2106c           ALE(i) = max(ale_wake(i),Ale_bl(i))
2107            ALE(i) = max(ale_wake(i),ale_bl_prescr)
2108con ne couple que alp
2109c           ALP(i) = alp_wake(i) + Alp_bl(i)
2110            ALP(i) = alp_wake(i) + alp_bl_prescr
2111          enddo
2112       else
2113         IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2114!         do i = 1,klon
2115!             ALE(i) = max(ale_wake(i),Ale_bl(i))
2116! avant        ALP(i) = alp_wake(i) + Alp_bl(i)
2117!             ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2118!         write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2119!         write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2120!         enddo
2121
2122!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2123! Modif FH 2010/04/27. Sans doute temporaire.
2124! Deux options pour le alp_offset : constant si >Â 0 ou proportionnel Ãa
2125! w si <0
2126!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2127       do i = 1,klon
2128          ALE(i) = max(ale_wake(i),Ale_bl(i))
2129          if (alp_offset>=0.) then
2130            ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2131          else
2132            ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
2133            if (alp(i)<0.) then
2134               print*,'ALP ',alp(i),alp_wake(i)
2135     s         ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
2136            endif
2137          endif
2138       enddo
2139!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2140
2141       endif
2142       do i=1,klon
2143          if (alp(i)>alp_max) then
2144             IF(prt_level>9)WRITE(lunout,*)                             &
2145     &       'WARNING SUPER ALP (seuil=',alp_max,
2146     ,       '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2147             alp(i)=alp_max
2148          endif
2149          if (ale(i)>ale_max) then
2150             IF(prt_level>9)WRITE(lunout,*)                             &
2151     &       'WARNING SUPER ALE (seuil=',ale_max,
2152     ,       '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2153             ale(i)=ale_max
2154          endif
2155       enddo
2156
2157cfin calcul ale et alp
2158c=================================================================================================
2159
2160
2161c sb, oct02:
2162c Schema de convection modularise et vectorise:
2163c (driver commun aux versions 3 et 4)
2164c
2165          IF (ok_cvl) THEN ! new driver for convectL
2166
2167          CALL concvl (iflag_con,iflag_clos,
2168     .        dtime,paprs,pplay,t_undi,q_undi,
2169     .        t_wake,q_wake,wake_s,
2170     .        u_seri,v_seri,tr_seri,nbtr,
2171     .        ALE,ALP,
2172     .        ema_work1,ema_work2,
2173     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2174     .        rain_con, snow_con, ibas_con, itop_con, sigd,
2175     .        ema_cbmf,upwd,dnwd,dnwd0,
2176     .        Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2177     .        pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2178     .        pmflxr,pmflxs,da,phi,mp,
2179     .        ftd,fqd,lalim_conv,wght_th)
2180
2181cIM begin
2182c       print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2183c    .dnwd0(1,1),ftd(1,1),fqd(1,1)
2184cIM end
2185cIM cf. FH
2186              clwcon0=qcondc
2187              pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2188
2189              do i = 1, klon
2190                if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2191              enddo
2192
2193          ELSE ! ok_cvl
2194
2195c MAF conema3 ne contient pas les traceurs
2196          CALL conema3 (dtime,
2197     .        paprs,pplay,t_seri,q_seri,
2198     .        u_seri,v_seri,tr_seri,ntra,
2199     .        ema_work1,ema_work2,
2200     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2201     .        rain_con, snow_con, ibas_con, itop_con,
2202     .        upwd,dnwd,dnwd0,bas,top,
2203     .        Ma,cape,tvp,rflag,
2204     .        pbase
2205     .        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2206     .        ,clwcon0)
2207
2208          ENDIF ! ok_cvl
2209
2210c
2211c Correction precip
2212          rain_con = rain_con * cvl_corr
2213          snow_con = snow_con * cvl_corr
2214c
2215
2216           IF (.NOT. ok_gust) THEN
2217           do i = 1, klon
2218            wd(i)=0.0
2219           enddo
2220           ENDIF
2221
2222c =================================================================== c
2223c Calcul des proprietes des nuages convectifs
2224c
2225
2226c   calcul des proprietes des nuages convectifs
2227             clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2228             call clouds_gno
2229     s       (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2230
2231c =================================================================== c
2232
2233          DO i = 1, klon
2234           itop_con(i) = min(max(itop_con(i),1),klev)
2235           ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2236          ENDDO
2237
2238          DO i = 1, klon
2239            ema_pcb(i)  = paprs(i,ibas_con(i))
2240          ENDDO
2241          DO i = 1, klon
2242! L'idicage de itop_con peut cacher un pb potentiel
2243! FH sous la dictee de JYG, CR
2244            ema_pct(i)  = paprs(i,itop_con(i)+1)
2245
2246            if (itop_con(i).gt.klev-3) then
2247              if(prt_level >= 9) then
2248                write(lunout,*)'La convection monte trop haut '
2249                write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2250              endif
2251            endif
2252          ENDDO     
2253      ELSE IF (iflag_con.eq.0) THEN
2254          write(lunout,*) 'On n appelle pas la convection'
2255          clwcon0=0.
2256          rnebcon0=0.
2257          d_t_con=0.
2258          d_q_con=0.
2259          d_u_con=0.
2260          d_v_con=0.
2261          rain_con=0.
2262          snow_con=0.
2263          bas=1
2264          top=1
2265      ELSE
2266          WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2267          CALL abort
2268      ENDIF
2269
2270c     CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2271c    .              d_u_con, d_v_con)
2272
2273!-----------------------------------------------------------------------------------------
2274! ajout des tendances de la diffusion turbulente
2275      CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2276!-----------------------------------------------------------------------------------------
2277
2278      if (mydebug) then
2279        call writefield_phy('u_seri',u_seri,llm)
2280        call writefield_phy('v_seri',v_seri,llm)
2281        call writefield_phy('t_seri',t_seri,llm)
2282        call writefield_phy('q_seri',q_seri,llm)
2283      endif
2284
2285cIM
2286      IF (ip_ebil_phy.ge.2) THEN
2287        ztit='after convect'
2288        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2289     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2290     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2291         call diagphy(airephy,ztit,ip_ebil_phy
2292     e      , zero_v, zero_v, zero_v, zero_v, zero_v
2293     e      , zero_v, rain_con, snow_con, ztsol
2294     e      , d_h_vcol, d_qt, d_ec
2295     s      , fs_bound, fq_bound )
2296      END IF
2297C
2298      IF (check) THEN
2299          za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2300          WRITE(lunout,*)"aprescon=", za
2301          zx_t = 0.0
2302          za = 0.0
2303          DO i = 1, klon
2304            za = za + airephy(i)/REAL(klon)
2305            zx_t = zx_t + (rain_con(i)+
2306     .                   snow_con(i))*airephy(i)/REAL(klon)
2307          ENDDO
2308          zx_t = zx_t/za*dtime
2309          WRITE(lunout,*)"Precip=", zx_t
2310      ENDIF
2311      IF (zx_ajustq) THEN
2312          DO i = 1, klon
2313            z_apres(i) = 0.0
2314          ENDDO
2315          DO k = 1, klev
2316            DO i = 1, klon
2317              z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2318     .            *(paprs(i,k)-paprs(i,k+1))/RG
2319            ENDDO
2320          ENDDO
2321          DO i = 1, klon
2322            z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2323     .          /z_apres(i)
2324          ENDDO
2325          DO k = 1, klev
2326            DO i = 1, klon
2327              IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2328     .            z_factor(i).LT.(1.0-1.0E-08)) THEN
2329                  q_seri(i,k) = q_seri(i,k) * z_factor(i)
2330              ENDIF
2331            ENDDO
2332          ENDDO
2333      ENDIF
2334      zx_ajustq=.FALSE.
2335
2336c
2337c=============================================================================
2338cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2339cpour la couche limite diffuse pour l instant
2340c
2341      if (iflag_wake.eq.1) then
2342      DO k=1,klev
2343        DO i=1,klon
2344          dt_dwn(i,k)  = ftd(i,k)
2345          wdt_PBL(i,k) = 0.
2346          dq_dwn(i,k)  = fqd(i,k)
2347          wdq_PBL(i,k) = 0.
2348          M_dwn(i,k)   = dnwd0(i,k)
2349          M_up(i,k)    = upwd(i,k)
2350          dt_a(i,k)    = d_t_con(i,k)/dtime - ftd(i,k)
2351          udt_PBL(i,k) = 0.
2352          dq_a(i,k)    = d_q_con(i,k)/dtime - fqd(i,k)
2353          udq_PBL(i,k) = 0.
2354        ENDDO
2355      ENDDO
2356c
2357ccalcul caracteristiques de la poche froide
2358      call calWAKE (paprs,pplay,dtime
2359     :               ,t_seri,q_seri,omega
2360     :               ,dt_dwn,dq_dwn,M_dwn,M_up
2361     :               ,dt_a,dq_a,sigd
2362     :               ,wdt_PBL,wdq_PBL
2363     :               ,udt_PBL,udq_PBL
2364     o               ,wake_deltat,wake_deltaq,wake_dth
2365     o               ,wake_h,wake_s,wake_dens
2366     o               ,wake_pe,wake_fip,wake_gfl
2367     o               ,dt_wake,dq_wake
2368     o               ,wake_k, t_undi,q_undi
2369     o               ,wake_omgbdth,wake_dp_omgb
2370     o               ,wake_dtKE,wake_dqKE
2371     o               ,wake_dtPBL,wake_dqPBL
2372     o               ,wake_omg,wake_dp_deltomg
2373     o               ,wake_spread,wake_Cstar,wake_d_deltat_gw
2374     o               ,wake_ddeltat,wake_ddeltaq)
2375c
2376!-----------------------------------------------------------------------------------------
2377! ajout des tendances des poches froides
2378! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2379! coherent avec les autres d_t_...
2380      d_t_wake(:,:)=dt_wake(:,:)*dtime
2381      d_q_wake(:,:)=dq_wake(:,:)*dtime
2382      CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2383!-----------------------------------------------------------------------------------------
2384
2385      endif
2386c
2387c===================================================================
2388cJYG
2389      IF (ip_ebil_phy.ge.2) THEN
2390        ztit='after wake'
2391        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2392     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2393     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2394        call diagphy(airephy,ztit,ip_ebil_phy
2395     e      , zero_v, zero_v, zero_v, zero_v, zero_v
2396     e      , zero_v, zero_v, zero_v, ztsol
2397     e      , d_h_vcol, d_qt, d_ec
2398     s      , fs_bound, fq_bound )
2399      END IF
2400
2401c      print*,'apres callwake iflag_cldcon=', iflag_cldcon
2402c
2403c===================================================================
2404c Convection seche (thermiques ou ajustement)
2405c===================================================================
2406c
2407       call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2408     s ,seuil_inversion,weak_inversion,dthmin)
2409
2410
2411
2412      d_t_ajsb(:,:)=0.
2413      d_q_ajsb(:,:)=0.
2414      d_t_ajs(:,:)=0.
2415      d_u_ajs(:,:)=0.
2416      d_v_ajs(:,:)=0.
2417      d_q_ajs(:,:)=0.
2418      clwcon0th(:,:)=0.
2419c
2420c      fm_therm(:,:)=0.
2421c      entr_therm(:,:)=0.
2422c      detr_therm(:,:)=0.
2423c
2424      IF(prt_level>9)WRITE(lunout,*)
2425     .    'AVANT LA CONVECTION SECHE , iflag_thermals='
2426     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
2427      if(iflag_thermals.lt.0) then
2428c  Rien
2429c  ====
2430         IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2431
2432
2433      else
2434
2435c  Thermiques
2436c  ==========
2437         IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2438     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
2439
2440
2441         if (iflag_thermals.gt.1) then
2442         call calltherm(pdtphys
2443     s      ,pplay,paprs,pphi,weak_inversion
2444     s      ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2445     s      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2446     s      ,fm_therm,entr_therm,detr_therm
2447     s      ,zqasc,clwcon0th,lmax_th,ratqscth
2448     s      ,ratqsdiff,zqsatth
2449con rajoute ale et alp, et les caracteristiques de la couche alim
2450     s      ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca
2451     s      ,ztv,zpspsk,ztla,zthl)
2452
2453! ----------------------------------------------------------------------
2454! Transport de la TKE par les panaches thermiques.
2455! FH : 2010/02/01
2456      if (iflag_pbl.eq.10) then
2457      call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
2458     s           rg,paprs,pbl_tke)
2459      endif
2460! ----------------------------------------------------------------------
2461
2462         endif
2463
2464
2465
2466c  Ajustement sec
2467c  ==============
2468
2469! Dans le cas où on active les thermiques, on fait partir l'ajustement
2470! a partir du sommet des thermiques.
2471! Dans le cas contraire, on demarre au niveau 1.
2472
2473         if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2474
2475         if(iflag_thermals.eq.0) then
2476            IF(prt_level>9)WRITE(lunout,*)'ajsec'
2477            limbas(:)=1
2478         else
2479            limbas(:)=lmax_th(:)
2480         endif
2481 
2482! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2483! pour des test de convergence numerique.
2484! Le nouveau ajsec est a priori mieux, meme pour le cas
2485! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2486! non nulles numeriquement pour des mailles non concernees.
2487
2488         if (iflag_thermals.eq.0) then
2489            CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2490     s      , d_t_ajsb, d_q_ajsb)
2491         else
2492            CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2493     s      , d_t_ajsb, d_q_ajsb)
2494         endif
2495
2496!-----------------------------------------------------------------------------------------
2497! ajout des tendances de l'ajustement sec ou des thermiques
2498      CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2499         d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2500         d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2501
2502!-----------------------------------------------------------------------------------------
2503
2504         endif
2505
2506      endif
2507c
2508c===================================================================
2509cIM
2510      IF (ip_ebil_phy.ge.2) THEN
2511        ztit='after dry_adjust'
2512        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2513     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2514     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2515        call diagphy(airephy,ztit,ip_ebil_phy
2516     e      , zero_v, zero_v, zero_v, zero_v, zero_v
2517     e      , zero_v, zero_v, zero_v, ztsol
2518     e      , d_h_vcol, d_qt, d_ec
2519     s      , fs_bound, fq_bound )
2520      END IF
2521
2522
2523c-------------------------------------------------------------------------
2524c  Caclul des ratqs
2525c-------------------------------------------------------------------------
2526
2527c      print*,'calcul des ratqs'
2528c   ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2529c   ----------------
2530c   on ecrase le tableau ratqsc calcule par clouds_gno
2531      if (iflag_cldcon.eq.1) then
2532         do k=1,klev
2533         do i=1,klon
2534            if(ptconv(i,k)) then
2535              ratqsc(i,k)=ratqsbas
2536     s        +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2537            else
2538               ratqsc(i,k)=0.
2539            endif
2540         enddo
2541         enddo
2542
2543c-----------------------------------------------------------------------
2544c  par nversion de la fonction log normale
2545c-----------------------------------------------------------------------
2546      else if (iflag_cldcon.eq.4) then
2547         ptconvth(:,:)=.false.
2548         ratqsc(:,:)=0.
2549         if(prt_level.ge.9) print*,'avant clouds_gno thermique'
2550         call clouds_gno
2551     s   (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2552         if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
2553       
2554       endif
2555
2556c   ratqs stables
2557c   -------------
2558
2559      if (iflag_ratqs.eq.0) then
2560
2561! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2562         do k=1,klev
2563            do i=1, klon
2564               ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2565     s         min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2566            enddo
2567         enddo
2568
2569! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2570! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2571! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2572! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2573! Il s'agit de differents tests dans la phase de reglage du modele
2574! avec thermiques.
2575
2576      else if (iflag_ratqs.eq.1) then
2577
2578         do k=1,klev
2579            do i=1, klon
2580               if (pplay(i,k).ge.60000.) then
2581                  ratqss(i,k)=ratqsbas
2582               else if ((pplay(i,k).ge.30000.).and.
2583     s            (pplay(i,k).lt.60000.)) then
2584                  ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2585     s            (60000.-pplay(i,k))/(60000.-30000.)
2586               else
2587                  ratqss(i,k)=ratqshaut
2588               endif
2589            enddo
2590         enddo
2591
2592      else
2593
2594         do k=1,klev
2595            do i=1, klon
2596               if (pplay(i,k).ge.60000.) then
2597                  ratqss(i,k)=ratqsbas
2598     s            *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2599               else if ((pplay(i,k).ge.30000.).and.
2600     s             (pplay(i,k).lt.60000.)) then
2601                    ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2602     s              (60000.-pplay(i,k))/(60000.-30000.)
2603               else
2604                    ratqss(i,k)=ratqshaut
2605               endif
2606            enddo
2607         enddo
2608      endif
2609
2610
2611
2612
2613c  ratqs final
2614c  -----------
2615
2616      if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2617     s    .or.iflag_cldcon.eq.4) then
2618
2619! On ajoute une constante au ratqsc*2 pour tenir compte de
2620! fluctuations turbulentes de petite echelle
2621
2622         do k=1,klev
2623            do i=1,klon
2624               if ((fm_therm(i,k).gt.1.e-10)) then
2625                  ratqsc(i,k)=sqrt(ratqsc(i,k)**2+0.05**2)
2626               endif
2627            enddo
2628         enddo
2629
2630!   les ratqs sont une combinaison de ratqss et ratqsc
2631       if(prt_level.ge.9)
2632     $       write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
2633
2634         if (tau_ratqs>1.e-10) then
2635            facteur=exp(-pdtphys/tau_ratqs)
2636         else
2637            facteur=0.
2638         endif
2639         ratqs(:,:)=ratqsc(:,:)*(1.-facteur)+ratqs(:,:)*facteur
2640!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2641! FH 22/09/2009
2642! La ligne ci-dessous faisait osciller le modele et donnait une solution
2643! assymptotique bidon et dépendant fortement du pas de temps.
2644!        ratqs(:,:)=sqrt(ratqs(:,:)**2+ratqss(:,:)**2)
2645!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2646         ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
2647      else
2648!   on ne prend que le ratqs stable pour fisrtilp
2649         ratqs(:,:)=ratqss(:,:)
2650      endif
2651
2652
2653c
2654c Appeler le processus de condensation a grande echelle
2655c et le processus de precipitation
2656c-------------------------------------------------------------------------
2657      CALL fisrtilp(dtime,paprs,pplay,
2658     .           t_seri, q_seri,ptconv,ratqs,
2659     .           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2660     .           rain_lsc, snow_lsc,
2661     .           pfrac_impa, pfrac_nucl, pfrac_1nucl,
2662     .           frac_impa, frac_nucl,
2663     .           prfl, psfl, rhcl,
2664     .           zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon )
2665
2666      WHERE (rain_lsc < 0) rain_lsc = 0.
2667      WHERE (snow_lsc < 0) snow_lsc = 0.
2668!-----------------------------------------------------------------------------------------
2669! ajout des tendances de la diffusion turbulente
2670      CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
2671!-----------------------------------------------------------------------------------------
2672      DO k = 1, klev
2673      DO i = 1, klon
2674         cldfra(i,k) = rneb(i,k)
2675         IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2676      ENDDO
2677      ENDDO
2678      IF (check) THEN
2679         za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2680         WRITE(lunout,*)"apresilp=", za
2681         zx_t = 0.0
2682         za = 0.0
2683         DO i = 1, klon
2684            za = za + airephy(i)/REAL(klon)
2685            zx_t = zx_t + (rain_lsc(i)
2686     .                  + snow_lsc(i))*airephy(i)/REAL(klon)
2687        ENDDO
2688         zx_t = zx_t/za*dtime
2689         WRITE(lunout,*)"Precip=", zx_t
2690      ENDIF
2691cIM
2692      IF (ip_ebil_phy.ge.2) THEN
2693        ztit='after fisrt'
2694        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2695     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2696     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2697        call diagphy(airephy,ztit,ip_ebil_phy
2698     e      , zero_v, zero_v, zero_v, zero_v, zero_v
2699     e      , zero_v, rain_lsc, snow_lsc, ztsol
2700     e      , d_h_vcol, d_qt, d_ec
2701     s      , fs_bound, fq_bound )
2702      END IF
2703
2704      if (mydebug) then
2705        call writefield_phy('u_seri',u_seri,llm)
2706        call writefield_phy('v_seri',v_seri,llm)
2707        call writefield_phy('t_seri',t_seri,llm)
2708        call writefield_phy('q_seri',q_seri,llm)
2709      endif
2710
2711c
2712c-------------------------------------------------------------------
2713c  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2714c-------------------------------------------------------------------
2715
2716c 1. NUAGES CONVECTIFS
2717c
2718cIM cf FH
2719c     IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2720      IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2721       snow_tiedtke=0.
2722c     print*,'avant calcul de la pseudo precip '
2723c     print*,'iflag_cldcon',iflag_cldcon
2724       if (iflag_cldcon.eq.-1) then
2725          rain_tiedtke=rain_con
2726       else
2727c       print*,'calcul de la pseudo precip '
2728          rain_tiedtke=0.
2729c         print*,'calcul de la pseudo precip 0'
2730          do k=1,klev
2731          do i=1,klon
2732             if (d_q_con(i,k).lt.0.) then
2733                rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2734     s         *(paprs(i,k)-paprs(i,k+1))/rg
2735             endif
2736          enddo
2737          enddo
2738       endif
2739c
2740c     call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2741c
2742
2743c Nuages diagnostiques pour Tiedtke
2744      CALL diagcld1(paprs,pplay,
2745cIM cf FH  .             rain_con,snow_con,ibas_con,itop_con,
2746     .             rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2747     .             diafra,dialiq)
2748      DO k = 1, klev
2749      DO i = 1, klon
2750      IF (diafra(i,k).GT.cldfra(i,k)) THEN
2751         cldliq(i,k) = dialiq(i,k)
2752         cldfra(i,k) = diafra(i,k)
2753      ENDIF
2754      ENDDO
2755      ENDDO
2756
2757      ELSE IF (iflag_cldcon.ge.3) THEN
2758c  On prend pour les nuages convectifs le max du calcul de la
2759c  convection et du calcul du pas de temps precedent diminue d'un facteur
2760c  facttemps
2761      facteur = pdtphys *facttemps
2762      do k=1,klev
2763         do i=1,klon
2764            rnebcon(i,k)=rnebcon(i,k)*facteur
2765            if (rnebcon0(i,k)*clwcon0(i,k).gt.rnebcon(i,k)*clwcon(i,k))
2766     s      then
2767                rnebcon(i,k)=rnebcon0(i,k)
2768                clwcon(i,k)=clwcon0(i,k)
2769            endif
2770         enddo
2771      enddo
2772
2773c
2774cjq - introduce the aerosol direct and first indirect radiative forcings
2775cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2776      IF (ok_ade.OR.ok_aie) THEN
2777         IF (.NOT. aerosol_couple)
2778     &        CALL readaerosol_optic(
2779     &        debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref,
2780     &        pdtphys, pplay, paprs, t_seri, rhcl, presnivs,
2781     &        mass_solu_aero, mass_solu_aero_pi,
2782     &        tau_aero, piz_aero, cg_aero,
2783     &        tausum_aero, tau3d_aero)
2784      ELSE
2785cIM 170310 BEG
2786         tausum_aero(:,:,:) = 0.
2787cIM 170310 END
2788         tau_aero(:,:,:,:) = 0.
2789         piz_aero(:,:,:,:) = 0.
2790         cg_aero(:,:,:,:)  = 0.
2791      ENDIF
2792
2793cIM calcul nuages par le simulateur ISCCP
2794c
2795#ifdef histISCCP
2796      IF (ok_isccp) THEN
2797c
2798cIM lecture invtau, tautab des fichiers formattes
2799c
2800      IF (debut) THEN
2801c$OMP MASTER
2802c
2803      open(99,file='tautab.formatted', FORM='FORMATTED')
2804      read(99,'(f30.20)') tautab_omp
2805      close(99)
2806c
2807      open(99,file='invtau.formatted',form='FORMATTED')
2808      read(99,'(i10)') invtau_omp
2809
2810c     print*,'calcul_simulISCCP invtau_omp',invtau_omp
2811c     write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
2812
2813      close(99)
2814c$OMP END MASTER
2815c$OMP BARRIER
2816      tautab=tautab_omp
2817      invtau=invtau_omp
2818c
2819      ENDIF !debut
2820c
2821cIM appel simulateur toutes les  NINT(freq_ISCCP/dtime) heures
2822       IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
2823#include "calcul_simulISCCP.h"
2824       ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
2825      ENDIF !ok_isccp
2826#endif
2827
2828c   On prend la somme des fractions nuageuses et des contenus en eau
2829
2830      if (iflag_cldcon>=5) then
2831
2832! Si on est sur un point touche par la convection profonde et pas
2833! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
2834! de la convection profonde.
2835
2836         print*,'TEST SCHEMA DE NUAGES '
2837         do k=1,klev
2838            do i=1,klon
2839               if (ptconv(i,k).and. .not. ptconvth(i,k)) then
2840                   cldfra(i,k)=rnebcon(i,k)
2841                   cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
2842               endif
2843            enddo
2844         enddo
2845      else
2846! Ancienne version
2847      cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2848      cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2849      endif
2850
2851      ENDIF
2852c
2853c 2. NUAGES STARTIFORMES
2854c
2855      IF (ok_stratus) THEN
2856      CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2857      DO k = 1, klev
2858      DO i = 1, klon
2859      IF (diafra(i,k).GT.cldfra(i,k)) THEN
2860         cldliq(i,k) = dialiq(i,k)
2861         cldfra(i,k) = diafra(i,k)
2862      ENDIF
2863      ENDDO
2864      ENDDO
2865      ENDIF
2866c
2867c Precipitation totale
2868c
2869      DO i = 1, klon
2870         rain_fall(i) = rain_con(i) + rain_lsc(i)
2871         snow_fall(i) = snow_con(i) + snow_lsc(i)
2872      ENDDO
2873cIM
2874      IF (ip_ebil_phy.ge.2) THEN
2875        ztit="after diagcld"
2876        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2877     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2878     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2879        call diagphy(airephy,ztit,ip_ebil_phy
2880     e      , zero_v, zero_v, zero_v, zero_v, zero_v
2881     e      , zero_v, zero_v, zero_v, ztsol
2882     e      , d_h_vcol, d_qt, d_ec
2883     s      , fs_bound, fq_bound )
2884      END IF
2885c
2886c Calculer l'humidite relative pour diagnostique
2887c
2888      DO k = 1, klev
2889      DO i = 1, klon
2890         zx_t = t_seri(i,k)
2891         IF (thermcep) THEN
2892            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2893            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2894            zx_qs  = MIN(0.5,zx_qs)
2895            zcor   = 1./(1.-retv*zx_qs)
2896            zx_qs  = zx_qs*zcor
2897         ELSE
2898           IF (zx_t.LT.t_coup) THEN
2899              zx_qs = qsats(zx_t)/pplay(i,k)
2900           ELSE
2901              zx_qs = qsatl(zx_t)/pplay(i,k)
2902           ENDIF
2903         ENDIF
2904         zx_rh(i,k) = q_seri(i,k)/zx_qs
2905         zqsat(i,k)=zx_qs
2906      ENDDO
2907      ENDDO
2908
2909cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2910c   equivalente a 2m (tpote) pour diagnostique
2911c
2912      DO i = 1, klon
2913       tpot(i)=zt2m(i)*(100000./paprs(i,1))**RKAPPA
2914       IF (thermcep) THEN
2915        IF(zt2m(i).LT.RTT) then
2916         Lheat=RLSTT
2917        ELSE
2918         Lheat=RLVTT
2919        ENDIF
2920       ELSE
2921        IF (zt2m(i).LT.RTT) THEN
2922         Lheat=RLSTT
2923        ELSE
2924         Lheat=RLVTT
2925        ENDIF
2926       ENDIF
2927       tpote(i) = tpot(i)*     
2928     . EXP((Lheat *qsat2m(i))/(RCPD*zt2m(i)))
2929      ENDDO
2930
2931      IF (config_inca /= 'none') THEN
2932#ifdef INCA
2933         CALL VTe(VTphysiq)
2934         CALL VTb(VTinca)
2935         calday = REAL(days_elapsed + 1) + jH_cur
2936
2937         call chemtime(itap+itau_phy-1, date0, dtime)
2938         IF (config_inca == 'aero') THEN
2939            CALL AEROSOL_METEO_CALC(
2940     $           calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2941     $           prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
2942         END IF
2943
2944         zxsnow_dummy(:) = 0.0
2945
2946         CALL chemhook_begin (calday,
2947     $                          days_elapsed+1,
2948     $                          jH_cur,
2949     $                          pctsrf(1,1),
2950     $                          rlat,
2951     $                          rlon,
2952     $                          airephy,
2953     $                          paprs,
2954     $                          pplay,
2955     $                          coefh,
2956     $                          pphi,
2957     $                          t_seri,
2958     $                          u,
2959     $                          v,
2960     $                          wo(:, :, 1),
2961     $                          q_seri,
2962     $                          zxtsol,
2963     $                          zxsnow_dummy,
2964     $                          solsw,
2965     $                          albsol1,
2966     $                          rain_fall,
2967     $                          snow_fall,
2968     $                          itop_con,
2969     $                          ibas_con,
2970     $                          cldfra,
2971     $                          iim,
2972     $                          jjm,
2973     $                          tr_seri,
2974     $                          ftsol,
2975     $                          paprs,
2976     $                          cdragh,
2977     $                          cdragm,
2978     $                          pctsrf,
2979     $                          pdtphys,
2980     $                            itap)
2981
2982         CALL VTe(VTinca)
2983         CALL VTb(VTphysiq)
2984#endif
2985      END IF !config_inca /= 'none'
2986c     
2987c Calculer les parametres optiques des nuages et quelques
2988c parametres pour diagnostiques:
2989c
2990
2991      IF (aerosol_couple) THEN
2992         mass_solu_aero(:,:)    = ccm(:,:,1)
2993         mass_solu_aero_pi(:,:) = ccm(:,:,2)
2994      END IF
2995
2996      if (ok_newmicro) then
2997      CALL newmicro (paprs, pplay,ok_newmicro,
2998     .            t_seri, cldliq, cldfra, cldtau, cldemi,
2999     .            cldh, cldl, cldm, cldt, cldq,
3000     .            flwp, fiwp, flwc, fiwc,
3001     e            ok_aie,
3002     e            mass_solu_aero, mass_solu_aero_pi,
3003     e            bl95_b0, bl95_b1,
3004     s            cldtaupi, re, fl, ref_liq, ref_ice)
3005      else
3006      CALL nuage (paprs, pplay,
3007     .            t_seri, cldliq, cldfra, cldtau, cldemi,
3008     .            cldh, cldl, cldm, cldt, cldq,
3009     e            ok_aie,
3010     e            mass_solu_aero, mass_solu_aero_pi,
3011     e            bl95_b0, bl95_b1,
3012     s            cldtaupi, re, fl)
3013     
3014      endif
3015c
3016c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3017c
3018      IF (MOD(itaprad,radpas).EQ.0) THEN
3019
3020      DO i = 1, klon
3021         albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
3022     .             + falb1(i,is_lic) * pctsrf(i,is_lic)
3023     .             + falb1(i,is_ter) * pctsrf(i,is_ter)
3024     .             + falb1(i,is_sic) * pctsrf(i,is_sic)
3025         albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
3026     .               + falb2(i,is_lic) * pctsrf(i,is_lic)
3027     .               + falb2(i,is_ter) * pctsrf(i,is_ter)
3028     .               + falb2(i,is_sic) * pctsrf(i,is_sic)
3029      ENDDO
3030
3031      if (mydebug) then
3032        call writefield_phy('u_seri',u_seri,llm)
3033        call writefield_phy('v_seri',v_seri,llm)
3034        call writefield_phy('t_seri',t_seri,llm)
3035       call writefield_phy('q_seri',q_seri,llm)
3036      endif
3037     
3038      IF (aerosol_couple) THEN
3039#ifdef INCA
3040         CALL radlwsw_inca
3041     e        (kdlon,kflev,dist, rmu0, fract, solaire,
3042     e        paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
3043     e        wo(:, :, 1),
3044     e        cldfra, cldemi, cldtau,
3045     s        heat,heat0,cool,cool0,radsol,albpla,
3046     s        topsw,toplw,solsw,sollw,
3047     s        sollwdown,
3048     s        topsw0,toplw0,solsw0,sollw0,
3049     s        lwdn0, lwdn, lwup0, lwup,
3050     s        swdn0, swdn, swup0, swup,
3051     e        ok_ade, ok_aie,
3052     e        tau_aero, piz_aero, cg_aero,
3053     s        topswad_aero, solswad_aero,
3054     s        topswad0_aero, solswad0_aero,
3055     s        topsw_aero, topsw0_aero,
3056     s        solsw_aero, solsw0_aero,
3057     e        cldtaupi,
3058     s        topswai_aero, solswai_aero)
3059           
3060#endif
3061      ELSE
3062
3063         CALL radlwsw
3064     e        (dist, rmu0, fract,
3065     e        paprs, pplay,zxtsol,albsol1, albsol2,
3066     e        t_seri,q_seri,wo,
3067     e        cldfra, cldemi, cldtau,
3068     e        ok_ade, ok_aie,
3069     e        tau_aero, piz_aero, cg_aero,
3070     e        cldtaupi,new_aod,
3071     e        zqsat, flwc, fiwc,
3072     s        heat,heat0,cool,cool0,radsol,albpla,
3073     s        topsw,toplw,solsw,sollw,
3074     s        sollwdown,
3075     s        topsw0,toplw0,solsw0,sollw0,
3076     s        lwdn0, lwdn, lwup0, lwup,
3077     s        swdn0, swdn, swup0, swup,
3078     s        topswad_aero, solswad_aero,
3079     s        topswai_aero, solswai_aero,
3080     o        topswad0_aero, solswad0_aero,
3081     o        topsw_aero, topsw0_aero,
3082     o        solsw_aero, solsw0_aero,
3083     o        topswcf_aero, solswcf_aero)
3084         
3085
3086      ENDIF ! aerosol_couple
3087      itaprad = 0
3088      ENDIF ! MOD(itaprad,radpas)
3089      itaprad = itaprad + 1
3090
3091      IF (iflag_radia.eq.0) THEN
3092         IF (prt_level.ge.9) THEN
3093            PRINT *,'--------------------------------------------------'
3094            PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3095            PRINT *,'>>>>           heat et cool mis a zero '
3096            PRINT *,'--------------------------------------------------'
3097         END IF
3098         heat=0.
3099         cool=0.
3100      END IF
3101
3102c
3103c Ajouter la tendance des rayonnements (tous les pas)
3104c
3105      DO k = 1, klev
3106      DO i = 1, klon
3107         t_seri(i,k) = t_seri(i,k)
3108     .               + (heat(i,k)-cool(i,k)) * dtime/RDAY
3109      ENDDO
3110      ENDDO
3111c
3112      if (mydebug) then
3113        call writefield_phy('u_seri',u_seri,llm)
3114        call writefield_phy('v_seri',v_seri,llm)
3115        call writefield_phy('t_seri',t_seri,llm)
3116        call writefield_phy('q_seri',q_seri,llm)
3117      endif
3118 
3119cIM
3120      IF (ip_ebil_phy.ge.2) THEN
3121        ztit='after rad'
3122        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3123     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3124     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3125        call diagphy(airephy,ztit,ip_ebil_phy
3126     e      , topsw, toplw, solsw, sollw, zero_v
3127     e      , zero_v, zero_v, zero_v, ztsol
3128     e      , d_h_vcol, d_qt, d_ec
3129     s      , fs_bound, fq_bound )
3130      END IF
3131c
3132c
3133c Calculer l'hydrologie de la surface
3134c
3135c      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3136c     .            agesno, ftsol,fqsurf,fsnow, ruis)
3137c
3138
3139c
3140c Calculer le bilan du sol et la derive de temperature (couplage)
3141c
3142      DO i = 1, klon
3143c         bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3144c a la demande de JLD
3145         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3146      ENDDO
3147c
3148cmoddeblott(jan95)
3149c Appeler le programme de parametrisation de l'orographie
3150c a l'echelle sous-maille:
3151c
3152      IF (ok_orodr) THEN
3153c
3154c  selection des points pour lesquels le shema est actif:
3155        igwd=0
3156        DO i=1,klon
3157        itest(i)=0
3158c        IF ((zstd(i).gt.10.0)) THEN
3159        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3160          itest(i)=1
3161          igwd=igwd+1
3162          idx(igwd)=i
3163        ENDIF
3164        ENDDO
3165c        igwdim=MAX(1,igwd)
3166c
3167        IF (ok_strato) THEN
3168       
3169          CALL drag_noro_strato(klon,klev,dtime,paprs,pplay,
3170     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
3171     e                   igwd,idx,itest,
3172     e                   t_seri, u_seri, v_seri,
3173     s                   zulow, zvlow, zustrdr, zvstrdr,
3174     s                   d_t_oro, d_u_oro, d_v_oro)
3175
3176       ELSE
3177        CALL drag_noro(klon,klev,dtime,paprs,pplay,
3178     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
3179     e                   igwd,idx,itest,
3180     e                   t_seri, u_seri, v_seri,
3181     s                   zulow, zvlow, zustrdr, zvstrdr,
3182     s                   d_t_oro, d_u_oro, d_v_oro)
3183       ENDIF
3184c
3185c  ajout des tendances
3186!-----------------------------------------------------------------------------------------
3187! ajout des tendances de la trainee de l'orographie
3188      CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
3189!-----------------------------------------------------------------------------------------
3190c
3191      ENDIF ! fin de test sur ok_orodr
3192c
3193      if (mydebug) then
3194        call writefield_phy('u_seri',u_seri,llm)
3195        call writefield_phy('v_seri',v_seri,llm)
3196        call writefield_phy('t_seri',t_seri,llm)
3197        call writefield_phy('q_seri',q_seri,llm)
3198      endif
3199     
3200      IF (ok_orolf) THEN
3201c
3202c  selection des points pour lesquels le shema est actif:
3203        igwd=0
3204        DO i=1,klon
3205        itest(i)=0
3206        IF ((zpic(i)-zmea(i)).GT.100.) THEN
3207          itest(i)=1
3208          igwd=igwd+1
3209          idx(igwd)=i
3210        ENDIF
3211        ENDDO
3212c        igwdim=MAX(1,igwd)
3213c
3214        IF (ok_strato) THEN
3215
3216          CALL lift_noro_strato(klon,klev,dtime,paprs,pplay,
3217     e                   rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval,
3218     e                   igwd,idx,itest,
3219     e                   t_seri, u_seri, v_seri,
3220     s                   zulow, zvlow, zustrli, zvstrli,
3221     s                   d_t_lif, d_u_lif, d_v_lif               )
3222       
3223        ELSE
3224          CALL lift_noro(klon,klev,dtime,paprs,pplay,
3225     e                   rlat,zmea,zstd,zpic,
3226     e                   itest,
3227     e                   t_seri, u_seri, v_seri,
3228     s                   zulow, zvlow, zustrli, zvstrli,
3229     s                   d_t_lif, d_u_lif, d_v_lif)
3230       ENDIF
3231c   
3232!-----------------------------------------------------------------------------------------
3233! ajout des tendances de la portance de l'orographie
3234      CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3235!-----------------------------------------------------------------------------------------
3236c
3237      ENDIF ! fin de test sur ok_orolf
3238C  HINES GWD PARAMETRIZATION
3239
3240       IF (ok_hines) then
3241
3242         CALL hines_gwd(klon,klev,dtime,paprs,pplay,
3243     i                  rlat,t_seri,u_seri,v_seri,
3244     o                  zustrhi,zvstrhi,
3245     o                  d_t_hin, d_u_hin, d_v_hin)
3246c
3247c  ajout des tendances
3248        CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'lif')
3249
3250      ENDIF
3251c
3252
3253c
3254cIM cf. FLott BEG
3255C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3256
3257      if (mydebug) then
3258        call writefield_phy('u_seri',u_seri,llm)
3259        call writefield_phy('v_seri',v_seri,llm)
3260        call writefield_phy('t_seri',t_seri,llm)
3261        call writefield_phy('q_seri',q_seri,llm)
3262      endif
3263
3264      DO i = 1, klon
3265        zustrph(i)=0.
3266        zvstrph(i)=0.
3267      ENDDO
3268      DO k = 1, klev
3269      DO i = 1, klon
3270       zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3271     c            (paprs(i,k)-paprs(i,k+1))/rg
3272       zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3273     c            (paprs(i,k)-paprs(i,k+1))/rg
3274      ENDDO
3275      ENDDO
3276c
3277cIM calcul composantes axiales du moment angulaire et couple des montagnes
3278c
3279      IF (is_sequential) THEN
3280     
3281        CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
3282     C                 ra,rg,romega,
3283     C                 rlat,rlon,pphis,
3284     C                 zustrdr,zustrli,zustrph,
3285     C                 zvstrdr,zvstrli,zvstrph,
3286     C                 paprs,u,v,
3287     C                 aam, torsfc)
3288       ENDIF
3289cIM cf. FLott END
3290cIM
3291      IF (ip_ebil_phy.ge.2) THEN
3292        ztit='after orography'
3293        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3294     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3295     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3296         call diagphy(airephy,ztit,ip_ebil_phy
3297     e      , zero_v, zero_v, zero_v, zero_v, zero_v
3298     e      , zero_v, zero_v, zero_v, ztsol
3299     e      , d_h_vcol, d_qt, d_ec
3300     s      , fs_bound, fq_bound )
3301      END IF
3302c
3303c
3304!====================================================================
3305! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
3306!====================================================================
3307! Abderrahmane 24.08.09
3308
3309      IF (ok_cosp) THEN
3310! adeclarer
3311#ifdef CPP_COSP
3312       IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
3313
3314       print*,'freq_cosp',freq_cosp
3315          mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
3316!       print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
3317!     s        ref_liq,ref_ice
3318          call phys_cosp(itap,dtime,freq_cosp,
3319     $                   ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP,
3320     $                   ecrit_mth,ecrit_day,ecrit_hf,
3321     $                   klon,klev,rlon,rlat,presnivs,overlap,
3322     $                   ref_liq,ref_ice,
3323     $                   pctsrf(:,is_ter)+pctsrf(:,is_lic),
3324     $                   zu10m,zv10m,
3325     $                   zphi,paprs(:,1:klev),pplay,zxtsol,t_seri,
3326     $                   qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc,
3327     $                   prfl(:,1:klev),psfl(:,1:klev),
3328     $                   pmflxr(:,1:klev),pmflxs(:,1:klev),
3329     $                   mr_ozone,cldtau, cldemi)
3330
3331!     L          calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
3332!     L          cfaddbze,clcalipso2,dbze,cltlidarradar,
3333!     M          clMISR,
3334!     R          clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
3335!     I          tauisccp,albisccp,meantbisccp,meantbclrisccp)
3336
3337         ENDIF
3338
3339#endif
3340       ENDIF  !ok_cosp
3341!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3342cAA
3343cAA Installation de l'interface online-offline pour traceurs
3344cAA
3345c====================================================================
3346c   Calcul  des tendances traceurs
3347c====================================================================
3348C
3349
3350      call phytrac (
3351     I     itap,     days_elapsed+1,    jH_cur,   debut,
3352     I     lafin,    dtime,     u, v,     t,
3353     I     paprs,    pplay,     pmfu,     pmfd,
3354     I     pen_u,    pde_u,     pen_d,    pde_d,
3355     I     cdragh,   coefh,     fm_therm, entr_therm,
3356     I     u1,       v1,        ftsol,    pctsrf,
3357     I     rlat,     frac_impa, frac_nucl,rlon,
3358     I     presnivs, pphis,     pphi,     albsol1,
3359     I     qx(:,:,ivap),rhcl,   cldfra,   rneb,
3360     I     diafra,   cldliq,    itop_con, ibas_con,
3361     I     pmflxr,   pmflxs,    prfl,     psfl,
3362     I     da,       phi,       mp,       upwd,     
3363     I     dnwd,     aerosol_couple,      flxmass_w,
3364     I     tau_aero, piz_aero,  cg_aero,  ccm,
3365     I     rfname,
3366     O     tr_seri)
3367
3368      IF (offline) THEN
3369
3370       IF (prt_level.ge.9)
3371     $    print*,'Attention on met a 0 les thermiques pour phystoke'
3372         call phystokenc (
3373     I                   nlon,klev,pdtphys,rlon,rlat,
3374     I                   t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3375     I                   fm_therm,entr_therm,
3376     I                   cdragh,coefh,u1,v1,ftsol,pctsrf,
3377     I                   frac_impa, frac_nucl,
3378     I                   pphis,airephy,dtime,itap)
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.