Context Navigation

physiq.F @ 1255

Last change on this file since 1255 was 1249, checked in by yann meurdesoif, 15 years ago
Corrections de Bug divers - portage vers Titane (CCRT) - YM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 115.2 KB

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

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

Download in other formats:

Original Format