Context Navigation

physiq.F @ 3536

Last change on this file since 3536 was 1311, checked in by Laurent Fairhead, 15 years ago

Modifications to thermals for TKE transport

Modifications aux thermiques pour le transport de la TKE

pbl_surface_mode.F90 : ok_flux_surf=.false. seulement pour klon>1
physiq.F : option iflag_pbl=10 pour transporter la TKE avec les thermiques.
calltherm.F90 : passage de iflag_thermals_ed en argument pour thermcell_plume
thermcell_main.F90 : Appel a plusieurs version de thermcell_plume en option
thermcell_plume.F90 : plusieurs versions dans le meme fichier (temporaire)
thermcell_height.F90 : verrue pour les cas ou les thermiques montent tout

en haut

yamada4 : inclusion de la diffusion verticale en option iflag_pbl=9

+ variables anciennement common, puis save/allocatable, remises en local

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 118.0 KB

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

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

Download in other formats:

Original Format