Context Navigation

physiq.F @ 1369

Last change on this file since 1369 was 1369, checked in by idelkadi, 14 years ago
Corrections pour integrer de nouveaux champs de sorties
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 122.6 KB

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

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

Download in other formats:

Original Format