Context Navigation

physiq.F @ 1299

Last change on this file since 1299 was 1299, checked in by Laurent Fairhead, 14 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

tous les FLOAT() sont remplacés par des REAL()
tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

all FLOAT() instructions have been replaced by REAL() instructions
all STOP instructions in phylmd have been replaced by calls to abort_gcm
the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
$Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

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

File size: 117.6 KB

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

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

Download in other formats:

Original Format