Context Navigation

physiq.F @ 1154

Last change on this file since 1154 was 1154, checked in by lguez, 15 years ago
-- Added "NetCDF95" interface in "bibio". -- NetCDF95 uses the module "typesizes", which is part of NetCDF, so we exclude dependency on "typesizes" in "bld.cfg". -- Added "assert_eq" and "assert" procedures, which are in the public part of Numerical Recipes. -- Added some interpolation and regridding utilities in "bibio". -- Added the ability to read an ozone climatology from a NetCDF file. -- Commented out unused variables and code in "etat0_netcdf". -- Updated calls to NetCDF in "etat0_netcdf": from Fortran 77 interface to Fortran 90 interface. -- Removed useless "deallocate" at the end of "etat0_netcdf". -- Corrected some declarations not conforming to Fortran standard, such as "integer4", or obsolescent such as "character4". -- Replaced some calls to not-standard function "float" by calls to "real". -- On Brodie at IDRIS, the NetCDF library compiled with OpenMP should be used. Changed path in "arch-SX8_BRODIE.path". -- Added warning for incompatibility of debugging options and OpenMP parallelization in "makelmdz_fcm".
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 116.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format