Context Navigation

physiq.F90 @ 2009

Last change on this file since 2009 was 2009, checked in by Laurent Fairhead, 10 years ago

Rajout de la clé CPP_RRTM pour la compilation avec RRTM

Added CPP_RRTM key for RRTM compilation

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 124.2 KB

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

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

Download in other formats:

Original Format