source: LMDZ5/branches/LMDZ5-dev032011/libf/phylmd/physiq.F @ 3241

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

Few more fields need to be set to zero when radiation is not activated


Mise à zéro de champs supplémentaires quand le rayonnement est désactivé

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