Context Navigation

physiq.F90 @ 1924

Last change on this file since 1924 was 1924, checked in by idelkadi, 10 years ago

Rajout de nouveaux diagnostiques pour la comparaison entre les fractions nuageuses sorties derictement dans LMDZ et celles sorties via le simulateur ISCCP. Routines concernees : physiq.F90, phys_output_ctrlout_mod.F90, phys_local_var_mod.F90 et phys_output_write_mod.F90

Rajout de la fraction d'ensoleillement par jour aux arguments de l'interface avec Cosp (entree pour le simulateur ISCCP). Routine concernée : physiq.F90

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: 119.8 KB

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

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

Download in other formats:

Original Format