Context Navigation

physiq.F90 @ 2339

Last change on this file since 2339 was 2298, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes -r2237:2291 into testing branch
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: 139.0 KB

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

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

Download in other formats:

Original Format