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