Context Navigation

physiq_mod.F @ 3519

Last change on this file since 3519 was 3451, checked in by emillour, 14 months ago
Venus PCM: Add "Age of Air" scheme from Maureen Cohen. Activated with flag "ok_aoa=y" and requires that there is also an "aoa" tracer in traceur.def MC
File size: 73.9 KB

Line
1	!
2	! $Id: $
3	!
4	MODULE physiq_mod
5
6	IMPLICIT NONE
7
8	CONTAINS
9
10	SUBROUTINE physiq (nlon,nlev,nqmax,
11	. debut,lafin,rjourvrai,gmtime,pdtphys,
12	. paprs,pplay,ppk,pphi,pphis,presnivs,
13	. u,v,t,qx,
14	. flxmw,
15	. d_u, d_v, d_t, d_qx, d_ps)
16
17	c======================================================================
18	c
19	c Modifications pour la physique de Venus
20	c S. Lebonnois (LMD/CNRS) Septembre 2005
21	c
22	c ---------------------------------------------------------------------
23	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
24	c
25	c Objet: Moniteur general de la physique du modele
26	cAA Modifications quant aux traceurs :
27	cAA - uniformisation des parametrisations ds phytrac
28	cAA - stockage des moyennes des champs necessaires
29	cAA en mode traceur off-line
30	c modif ( P. Le Van , 12/10/98 )
31	c
32	c Arguments:
33	c
34	c nlon----input-I-nombre de points horizontaux
35	c nlev----input-I-nombre de couches verticales
36	c nqmax---input-I-nombre de traceurs
37	c debut---input-L-variable logique indiquant le premier passage
38	c lafin---input-L-variable logique indiquant le dernier passage
39	c rjour---input-R-numero du jour de l'experience
40	c gmtime--input-R-fraction de la journee (0 a 1)
41	c pdtphys-input-R-pas d'integration pour la physique (seconde)
42	c paprs---input-R-pression pour chaque inter-couche (en Pa)
43	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
44	c ppk ---input-R-fonction d'Exner au milieu de couche
45	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
46	c pphis---input-R-geopotentiel du sol
47	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
48	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
49	c v-------input-R-vitesse Y (de S a N) en m/s
50	c t-------input-R-temperature (K)
51	c qx------input-R-mass mixing ratio traceurs (kg/kg)
52	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
53	c flxmw---input-R-flux de masse vertical en kg/s
54	c
55	c d_u-----output-R-tendance physique de "u" (m/s/s)
56	c d_v-----output-R-tendance physique de "v" (m/s/s)
57	c d_t-----output-R-tendance physique de "t" (K/s)
58	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
59	c d_ps----output-R-tendance physique de la pression au sol
60	c======================================================================
61	USE ioipsl
62	! USE histcom ! not needed; histcom is included in ioipsl
63	use dimphy
64	USE geometry_mod,only: longitude, latitude, ! in radians
65	& longitude_deg,latitude_deg, ! in degrees
66	& cell_area,dx,dy
67	USE phys_state_var_mod ! Variables sauvegardees de la physique
68	USE cpdet_phy_mod, only: cpdet, t2tpot
69	USE chemparam_mod
70	USE age_of_air_mod, ONLY: ok_aoa, reinit_aoa, i_aoa, init_age
71	USE age_of_air_mod, ONLY: aoa_ini, age_of_air
72	USE conc
73	USE param_v4_h
74	USE compo_hedin83_mod2
75	use radlwsw_newtoncool_mod, only: radlwsw_newtoncool
76	! use ieee_arithmetic
77	use time_phylmdz_mod, only: annee_ref, day_ref, itau_phy
78	use mod_grid_phy_lmdz, only: nbp_lon
79	use infotrac_phy, only: iflag_trac, tname, ttext
80	use vertical_layers_mod, only: pseudoalt
81	use turb_mod, only : sens, turb_resolved
82	use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
83	use sed_and_prod_mad, only: aer_sedimentation, drop_sedimentation
84	use iono_h, only: temp_elect, temp_ion
85	#ifdef CPP_XIOS
86	use xios_output_mod, only: initialize_xios_output,
87	& update_xios_timestep,
88	& send_xios_field
89	use wxios, only: wxios_context_init, xios_context_finalize
90	#endif
91	#ifdef MESOSCALE
92	use comm_wrf
93	#else
94	use iophy
95	use write_field_phy
96	use mod_phys_lmdz_omp_data, ONLY: is_omp_master
97	USE mod_phys_lmdz_para, only : is_parallel,jj_nb,
98	& is_north_pole_phy,
99	& is_south_pole_phy,
100	& is_master
101	#endif
102	IMPLICIT none
103	c======================================================================
104	c CLEFS CPP POUR LES IO
105	c =====================
106	#ifndef MESOSCALE
107	c#define histhf
108	#define histday
109	#define histmth
110	#define histins
111	#endif
112	c======================================================================
113	#include "dimsoil.h"
114	#include "clesphys.h"
115	#include "iniprint.h"
116	#include "timerad.h"
117	#include "tabcontrol.h"
118	#include "nirdata.h"
119	#include "hedin.h"
120	c======================================================================
121	LOGICAL ok_journe ! sortir le fichier journalier
122	save ok_journe
123	c PARAMETER (ok_journe=.true.)
124	c
125	LOGICAL ok_mensuel ! sortir le fichier mensuel
126	save ok_mensuel
127	c PARAMETER (ok_mensuel=.true.)
128	c
129	LOGICAL ok_instan ! sortir le fichier instantane
130	save ok_instan
131	c PARAMETER (ok_instan=.true.)
132	c
133	c======================================================================
134	c
135	c Variables argument:
136	c
137	INTEGER nlon
138	INTEGER nlev
139	INTEGER nqmax
140	REAL rjourvrai
141	REAL gmtime
142	REAL pdtphys
143	LOGICAL debut, lafin
144	REAL paprs(klon,klev+1)
145	REAL pplay(klon,klev)
146	REAL pphi(klon,klev)
147	REAL pphis(klon)
148	REAL presnivs(klev)
149
150	! ADAPTATION GCM POUR CP(T)
151	REAL ppk(klon,klev)
152
153	REAL u(klon,klev)
154	REAL v(klon,klev)
155	REAL t(klon,klev)
156	REAL qx(klon,klev,nqmax)
157
158	REAL d_u_dyn(klon,klev)
159	REAL d_t_dyn(klon,klev)
160
161	REAL flxmw(klon,klev)
162
163	REAL d_u(klon,klev)
164	REAL d_v(klon,klev)
165	REAL d_t(klon,klev)
166	REAL d_qx(klon,klev,nqmax)
167	REAL d_ps(klon)
168
169	logical ok_hf
170	real ecrit_hf
171	integer nid_hf
172	save ok_hf, ecrit_hf, nid_hf
173
174	#ifdef histhf
175	data ok_hf,ecrit_hf/.true.,0.25/
176	#else
177	data ok_hf/.false./
178	#endif
179
180	c Variables propres a la physique
181
182	integer,save :: itap ! physics counter
183	REAL delp(klon,klev) ! epaisseur d'une couche
184	REAL omega(klon,klev) ! vitesse verticale en Pa/s (+ downward)
185	REAL vertwind(klon,klev) ! vitesse verticale en m/s (+ upward)
186
187	INTEGER igwd,idx(klon),itest(klon)
188
189	c Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro
190
191	REAL zulow(klon),zvlow(klon)
192	REAL zustrdr(klon), zvstrdr(klon)
193	REAL zustrli(klon), zvstrli(klon)
194	REAL zustrhi(klon), zvstrhi(klon)
195	REAL zublstrdr(klon), zvblstrdr(klon)
196	REAL znlow(klon), zeff(klon)
197	REAL zbl(klon), knu2(klon),kbreak(nlon)
198	REAL tau0(klon), ztau(klon,klev)
199
200	c Pour calcul GW drag oro et nonoro: CALCUL de N2:
201	real zdtlev(klon,klev),zdzlev(klon,klev)
202	real ztlev(klon,klev)
203	real zn2(klon,klev) ! BV^2 at plev
204
205	c Pour les bilans de moment angulaire,
206	integer bilansmc
207	c Pour le transport de ballons
208	integer ballons
209	c j'ai aussi besoin
210	c du stress de couche limite a la surface:
211
212	REAL zustrcl(klon),zvstrcl(klon)
213
214	c et du stress total c de la physique:
215
216	REAL zustrph(klon),zvstrph(klon)
217
218	c Variables locales:
219	c
220	REAL cdragh(klon) ! drag coefficient pour T and Q
221	REAL cdragm(klon) ! drag coefficient pour vent
222	c
223	cAA Pour TRACEURS
224	cAA
225	REAL,save,allocatable :: source(:,:)
226	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
227	REAL yu1(klon) ! vents dans la premiere couche U
228	REAL yv1(klon) ! vents dans la premiere couche V
229
230	REAL dsens(klon) ! derivee chaleur sensible
231	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
232	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
233	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
234	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
235	c
236	REAL Fsedim(klon,klev+1) ! Flux de sedimentation (kg.m-2)
237
238	c======================================================================
239	c
240	c Declaration des procedures appelees
241	c
242	EXTERNAL ajsec ! ajustement sec
243	EXTERNAL clmain ! couche limite
244	EXTERNAL clmain_ideal ! couche limite simple
245	EXTERNAL hgardfou ! verifier les temperatures
246	c EXTERNAL orbite ! calculer l'orbite
247	EXTERNAL phyetat0 ! lire l'etat initial de la physique
248	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
249	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
250	! EXTERNAL suphec ! initialiser certaines constantes
251	EXTERNAL transp ! transport total de l'eau et de l'energie
252	EXTERNAL printflag
253	EXTERNAL zenang
254	EXTERNAL diagetpq
255	EXTERNAL conf_phys
256	EXTERNAL diagphy
257	EXTERNAL mucorr
258	EXTERNAL nirco2abs
259	EXTERNAL nir_leedat
260	EXTERNAL nltecool
261	EXTERNAL nlte_tcool
262	EXTERNAL nlte_setup
263	EXTERNAL blendrad
264	EXTERNAL nlthermeq
265	EXTERNAL euvheat
266	EXTERNAL param_read_e107
267	EXTERNAL conduction
268	EXTERNAL molvis
269	EXTERNAL moldiff_red
270
271	c
272	c Variables locales
273	c
274	CXXX PB
275	REAL fluxt(klon,klev) ! flux turbulent de chaleur
276	REAL fluxu(klon,klev) ! flux turbulent de vitesse u
277	REAL fluxv(klon,klev) ! flux turbulent de vitesse v
278	c
279	REAL flux_dyn(klon,klev) ! flux de chaleur produit par la dynamique
280	REAL flux_ajs(klon,klev) ! flux de chaleur ajustement sec
281	REAL flux_ec(klon,klev) ! flux de chaleur Ec
282	c
283	REAL tmpout(klon,klev) ! [K/s]
284	c
285	REAL dist, rmu0(klon), fract(klon)
286	REAL zdtime, zlongi
287	c
288	INTEGER i, k, iq, ig, j, ll, ilon, ilat, ilev, isoil
289	c
290	REAL zphi(klon,klev)
291	REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
292	REAL tlaymean ! valeur temporaire pour calculer zzlay
293	real tsurf(klon)
294
295	c va avec nlte_tcool
296	INTEGER ierr_nlte
297	REAL varerr
298
299	! photochemistry
300
301	integer :: chempas
302	real :: zctime
303
304	! sedimentation
305
306	REAL :: m0_mode1(klev,2),m0_mode2(klev,2)
307	REAL :: m3_mode1(klev,3),m3_mode2 (klev,3)
308	REAL :: d_drop_sed(klev),d_ccn_sed(klev,2),d_liq_sed(klev,2)
309	REAL :: aer_flux(klev)
310	c
311	c Variables du changement
312	c
313	c ajs: ajustement sec
314	c vdf: couche limite (Vertical DiFfusion)
315	REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax)
316	REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
317	c
318	REAL d_ts(klon)
319	c
320	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
321	REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax)
322	c
323	CMOD LOTT: Tendances Orography Sous-maille
324	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
325	REAL d_t_oro(klon,klev)
326	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
327	REAL d_t_lif(klon,klev)
328	C Tendances Ondes de G non oro (runs strato).
329	REAL d_u_hin(klon,klev), d_v_hin(klon,klev)
330	REAL d_t_hin(klon,klev)
331
332	c Tendencies due to radiative scheme [K/s]
333	c d_t_rad,dtsw,dtlw,d_t_nirco2,d_t_nlte,d_t_euv
334	c are not computed at each physical timestep
335	c therefore, they are defined and saved in phys_state_var_mod
336
337	c Tendencies due to molecular viscosity and conduction
338	real d_t_conduc(klon,klev) ! [K/s]
339	real d_u_molvis(klon,klev) ! [m/s] /s
340	real d_v_molvis(klon,klev) ! [m/s] /s
341
342	c Tendencies due to molecular diffusion
343	real d_q_moldif(klon,klev,nqmax)
344
345	c Tendencies due to ambipolar ion diffusion
346	real d_q_iondif(klon,klev,nqmax)
347
348	c
349	c Variables liees a l'ecriture de la bande histoire physique
350	c
351	INTEGER ecrit_mth
352	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
353	c
354	INTEGER ecrit_day
355	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
356	c
357	INTEGER ecrit_ins
358	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
359	c
360	integer itau_w ! pas de temps ecriture = itap + itau_phy
361
362	c Variables locales pour effectuer les appels en serie
363	c
364	REAL t_seri(klon,klev)
365	REAL u_seri(klon,klev), v_seri(klon,klev)
366	c
367	REAL :: tr_seri(klon,klev,nqmax)
368	REAL :: tr_hedin(klon,klev,nqmax)
369	REAL :: d_tr(klon,klev,nqmax)
370	c pour sorties
371	REAL :: col_dens_tr(klon,nqmax)
372	REAL,allocatable,save :: prod_tr(:,:,:)
373	REAL,allocatable,save :: loss_tr(:,:,:)
374
375	c pour ioipsl
376	INTEGER nid_day, nid_mth, nid_ins
377	SAVE nid_day, nid_mth, nid_ins
378	INTEGER nhori, nvert, idayref
379	REAL zsto, zout, zsto1, zsto2, zero
380	parameter (zero=0.0e0)
381	real zjulian
382	save zjulian
383
384	CHARACTER*2 str2
385	character*20 modname
386	character*80 abort_message
387	logical ok_sync
388
389	character*30 nom_fichier
390	character*10 varname
391	character*40 vartitle
392	character*20 varunits
393	C Variables liees au bilan d'energie et d'enthalpi
394	REAL ztsol(klon)
395	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
396	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
397	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
398	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
399	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
400	REAL d_h_vcol_phy
401	REAL fs_bound, fq_bound
402	SAVE d_h_vcol_phy
403	REAL zero_v(klon),zero_v2(klon,klev)
404	CHARACTER*15 ztit
405	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
406	SAVE ip_ebil
407	DATA ip_ebil/2/
408	INTEGER if_ebil ! level for energy conserv. dignostics
409	SAVE if_ebil
410	c+jld ec_conser
411	REAL d_t_ec(klon,klev) ! tendance du a la conversion Ec -> E thermique
412	c-jld ec_conser
413
414	c ALBEDO VARIATIONS (VCD)
415	REAL factAlb
416	c TEST VENUS...
417	REAL mang(klon,klev) ! moment cinetique
418	REAL mangtot ! moment cinetique total
419
420	c cell_area for outputs in hist*
421	REAL cell_area_out(klon)
422	#ifdef MESOSCALE
423	REAL :: dt_dyn(klev)
424	#endif
425	c Declaration des constantes et des fonctions thermodynamiques
426	c
427	#include "YOMCST.h"
428
429	c======================================================================
430	c======================================================================
431	c INITIALISATIONS
432	c======================================================================
433
434	modname = 'physiq'
435	ok_sync=.TRUE.
436
437	bilansmc = 0
438	ballons = 0
439	! NE FONCTIONNENT PAS ENCORE EN PARALLELE !!!
440	#ifndef MESOSCALE
441	if (is_parallel) then
442	bilansmc = 0
443	ballons = 0
444	endif
445	#endif
446	IF (if_ebil.ge.1) THEN
447	DO i=1,klon
448	zero_v(i)=0.
449	END DO
450	DO i=1,klon
451	DO j=1,klev
452	zero_v2(i,j)=0.
453	END DO
454	END DO
455	END IF
456
457	c======================================================================
458	c DONE ONLY AT FIRST CALL
459	c========================
460	IF (debut) THEN
461	allocate(source(klon,nqmax))
462	allocate(prod_tr(klon,klev,nqmax))
463	allocate(loss_tr(klon,klev,nqmax))
464	allocate(no_emission(klon,klev))
465	allocate(o2_emission(klon,klev))
466
467	#ifdef CPP_XIOS
468	! Initialize XIOS context
469	write(,) "physiq: call wxios_context_init"
470	CALL wxios_context_init
471	#endif
472
473	! The call to suphec is now done in iniphysiq_mod (interface)
474	! CALL suphec ! initialiser constantes et parametres phys.
475
476	IF (if_ebil.ge.1) d_h_vcol_phy=0.
477	!
478	! load flag and parameter values from physiq.def
479	!
480	call conf_phys(ok_journe, ok_mensuel,
481	. ok_instan,
482	. if_ebil)
483
484	call phys_state_var_init(nqmax)
485	c
486	c Initialising Hedin model for upper atm
487	c (to be revised when coupled to chemistry) :
488	call conc_init
489
490	! initialise physics counter
491
492	itap = 0
493
494	#ifdef MESOSCALE
495	print*,'check pdtphys',pdtphys
496	PRINT*,'check phisfi ',pphis(1),pphis(klon)
497	PRINT*,'check geop',pphi(1,1),pphi(klon,klev)
498	PRINT*,'check radsol',radsol(1),radsol(klon)
499	print*,'check ppk',ppk(1,1),ppk(klon,klev)
500	print*,'check ftsoil',ftsoil(1,1),ftsoil(klon,nsoilmx)
501	print*,'check ftsol',ftsol(1),ftsol(klon)
502	print*, "check temp", t(1,1),t(klon,klev)
503	print*, "check pres",paprs(1,1),paprs(klon,klev),pplay(1,1),
504	. pplay(klon,klev)
505	print*, "check u", u(1,1),u(klon,klev)
506	print*, "check v", v(1,1),v(klon,klev)
507	print*,'check falbe',falbe(1),falbe(klon)
508	!nqtot=nqmax
509	!ALLOCATE(tname(nqtot))
510	!tname=noms
511	zmea=0.
512	zstd=0.
513	zsig=0.
514	zgam=0.
515	zthe=0.
516	dtime=pdtphys
517	#else
518	c
519	c Lecture startphy.nc :
520	c
521	CALL phyetat0 ("startphy.nc")
522	IF (.not.startphy_file) THEN
523	! Additionnal academic initializations
524	ftsol(:)=t(:,1) ! surface temperature as in first atm. layer
525	DO isoil=1, nsoilmx
526	! subsurface temperatures equal to surface temperature
527	ftsoil(:,isoil)=ftsol(:)
528	ENDDO
529	ENDIF
530	#endif
531
532	c dtime est defini dans tabcontrol.h et lu dans startphy
533	c pdtphys est calcule a partir des nouvelles conditions:
534	c Reinitialisation du pas de temps physique quand changement
535	IF (ABS(dtime-pdtphys).GT.0.001) THEN
536	WRITE(lunout,*) 'Pas physique a change',dtime,
537	. pdtphys
538	c abort_message='Pas physique n est pas correct '
539	c call abort_physic(modname,abort_message,1)
540	c----------------
541	c pour initialiser convenablement le time_counter, il faut tenir compte
542	c du changement de dtime en changeant itau_phy (point de depart)
543	itau_phy = NINT(itau_phy*dtime/pdtphys)
544	c----------------
545	dtime=pdtphys
546	ENDIF
547
548	radpas = NINT(RDAY/pdtphys/nbapp_rad)
549
550	CALL printflag( ok_journe,ok_instan )
551
552	#ifdef CPP_XIOS
553	write(,) "physiq: call initialize_xios_output"
554	call initialize_xios_output(rjourvrai,gmtime,pdtphys,RDAY,
555	& presnivs,pseudoalt)
556	#endif
557
558	c
559	c---------
560	c FLOTT
561	IF (ok_orodr) THEN
562	DO i=1,klon
563	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
564	ENDDO
565	CALL SUGWD(klon,klev,paprs,pplay)
566	DO i=1,klon
567	zuthe(i)=0.
568	zvthe(i)=0.
569	if(zstd(i).gt.10.)then
570	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
571	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
572	endif
573	ENDDO
574	ENDIF
575
576	if (bilansmc.eq.1) then
577	C OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
578	C DU BILAN DE MOMENT ANGULAIRE.
579	open(27,file='aaam_bud.out',form='formatted')
580	open(28,file='fields_2d.out',form='formatted')
581	write(,)'Ouverture de aaam_bud.out (FL Vous parle)'
582	write(,)'Ouverture de fields_2d.out (FL Vous parle)'
583	endif !bilansmc
584
585	c--------------SLEBONNOIS
586	C OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
587	C DES BALLONS
588	if (ballons.eq.1) then
589	open(30,file='ballons-lat.out',form='formatted')
590	open(31,file='ballons-lon.out',form='formatted')
591	open(32,file='ballons-u.out',form='formatted')
592	open(33,file='ballons-v.out',form='formatted')
593	open(34,file='ballons-alt.out',form='formatted')
594	write(,)'Ouverture des ballons*.out'
595	endif !ballons
596	c-------------
597
598	c---------
599	C TRACEURS
600	C source dans couche limite
601	source(:,:) = 0.0 ! pas de source, pour l'instant
602	c---------
603
604	c---------
605	c INITIALIZE THERMOSPHERIC PARAMETERS
606	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
607
608	if (callthermos) then
609	call fill_data_thermos
610	call allocate_param_thermos(klev)
611	call param_read_e107
612	endif
613
614	c Initialisation (recomputed in concentration2)
615	do ig=1,klon
616	do j=1,klev
617	rnew(ig,j)=RD
618	cpnew(ig,j)=cpdet(t(ig,j))
619	mmean(ig,j)=RMD
620	akknew(ig,j)=1.e-4
621	rho(ig,j)=pplay(ig,j)/(rnew(ig,j)*t(ig,j))
622	enddo
623
624	enddo
625
626	IF(callthermos.or.callnlte.or.callnirco2) THEN
627	call compo_hedin83_init2
628	ENDIF
629	if (callnlte.and.nltemodel.eq.2) call nlte_setup
630	if (callnirco2.and.nircorr.eq.1) call nir_leedat
631	c---------
632
633	c
634	c Verifications:
635	c
636	IF (nlon .NE. klon) THEN
637	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
638	. klon
639	abort_message='nlon et klon ne sont pas coherents'
640	call abort_physic(modname,abort_message,1)
641	ENDIF
642	IF (nlev .NE. klev) THEN
643	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
644	. klev
645	abort_message='nlev et klev ne sont pas coherents'
646	call abort_physic(modname,abort_message,1)
647	ENDIF
648	c
649	IF (dtimeREAL(radpas).GT.(RDAY0.25).AND.cycle_diurne)
650	$ THEN
651	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
652	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
653	abort_message='Nbre d appels au rayonnement insuffisant'
654	call abort_physic(modname,abort_message,1)
655	ENDIF
656	c
657	WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
658	. iflag_ajs
659	c
660	ecrit_mth = NINT(RDAY/dtime*ecriphy) ! tous les ecritphy jours
661	IF (ok_mensuel) THEN
662	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
663	. ecrit_mth
664	ENDIF
665
666	ecrit_day = NINT(RDAY/dtime *1.0) ! tous les jours
667	IF (ok_journe) THEN
668	WRITE(lunout,*)'La frequence de sortie journaliere est de ',
669	. ecrit_day
670	ENDIF
671
672	ecrit_ins = NINT(RDAY/dtime*ecriphy) ! Fraction de jour reglable
673	IF (ok_instan) THEN
674	WRITE(lunout,*)'La frequence de sortie instant. est de ',
675	. ecrit_ins
676	ENDIF
677
678	c Verification synchronize AMBIPOLAR DIFFUSION & CHEMISTRY
679
680	IF ((ok_iondiff) .and. (NINT(RDAY/dtime).ne.nbapp_chem)) THEN
681	WRITE(lunout,*)'nbapp_chem .NE. day_step'
682	WRITE(lunout,*)'nbapp_chem ', nbapp_chem
683	WRITE(lunout,*)'day_step ', NINT(RDAY/dtime)
684	WRITE(lunout,*)'nbapp_chem must be equal to day_step'
685	STOP
686	ENDIF
687
688	c Initialisation des sorties
689	c===========================
690
691	#ifdef CPP_IOIPSL
692
693	#ifdef histhf
694	#include "ini_histhf.h"
695	#endif
696
697	#ifdef histday
698	#include "ini_histday.h"
699	#endif
700
701	#ifdef histmth
702	#include "ini_histmth.h"
703	#endif
704
705	#ifdef histins
706	#include "ini_histins.h"
707	#endif
708
709	#endif
710
711	c
712	c Initialiser les valeurs de u pour calculs tendances
713	c (pour T, c'est fait dans phyetat0)
714	c
715	DO k = 1, klev
716	DO i = 1, klon
717	u_ancien(i,k) = u(i,k)
718	ENDDO
719	ENDDO
720
721	! initialisation of microphysical and chemical parameters
722
723	if (ok_chem .and. .not. ok_cloud) then
724	print*,"chemistry requires clouds"
725	print*,"ok_cloud must be .true."
726	stop
727	end if
728
729	if (.not. ok_chem .and. ok_cloud .and. (cl_scheme == 1)) then
730	print*,"cl_scheme = 1 requires chemistry"
731	print*,"ok_chem must be .true."
732	stop
733	end if
734
735	! number of microphysical tracers
736
737	nmicro = 0
738	if (ok_cloud .and. (cl_scheme == 1)) nmicro = 2
739	if (ok_cloud .and. (cl_scheme == 2)) nmicro = 12
740
741	! initialise chemical parameters. includes the indexation of microphys tracers
742
743	if (ok_chem .or. cl_scheme == 2) then
744	call chemparam_ini()
745	end if
746
747	if ((iflag_trac.eq.1) .and. ok_aoa) then
748	write(lunout,*) 'Initialising age of air subroutine'
749	allocate(init_age(klon,klev))
750	call aoa_ini(ok_chem, tr_scheme) ! get index of age of air tracer in the tracer array
751	if (reinit_aoa) then
752	age(:,:) = 0. ! Set initial value of age of air to 0 everywhere
753	tr_seri(:,:,i_aoa) = 1e-30 ! Set initial value of tracer to tiny everywhere
754	endif ! reinitialisation loop
755	init_age(:,:) = age(:,:) ! save initial value of age, either read in from restartphy or set to 0
756	endif ! age of air initialisation
757
758	! initialise cloud parameters (for cl_scheme = 1)
759
760	if (ok_cloud .and. (cl_scheme == 1)) then
761	call cloud_ini(nlon,nlev,nb_mode)
762	end if
763
764	! initialise mmean
765
766	if (callthermos) then
767	call concentrations2(pplay,t,qx,nqmax)
768	end if
769
770	c========================
771	ENDIF ! debut
772	c========================
773
774	c======================================================================
775	! ------------------------------------------------------
776	! Initializations done at every physical timestep:
777	! ------------------------------------------------------
778
779	c Mettre a zero des variables de sortie (pour securite)
780	c
781	DO i = 1, klon
782	d_ps(i) = 0.0
783	ENDDO
784	DO k = 1, klev
785	DO i = 1, klon
786	d_t(i,k) = 0.0
787	d_u(i,k) = 0.0
788	d_v(i,k) = 0.0
789	ENDDO
790	ENDDO
791	DO iq = 1, nqmax
792	DO k = 1, klev
793	DO i = 1, klon
794	d_qx(i,k,iq) = 0.0
795	ENDDO
796	ENDDO
797	ENDDO
798	c
799	c Ne pas affecter les valeurs entrees de u, v, h, et q
800	c
801	DO k = 1, klev
802	DO i = 1, klon
803	t_seri(i,k) = t(i,k)
804	u_seri(i,k) = u(i,k)
805	v_seri(i,k) = v(i,k)
806	ENDDO
807	ENDDO
808	DO iq = 1, nqmax
809	DO k = 1, klev
810	DO i = 1, klon
811	tr_seri(i,k,iq) = qx(i,k,iq)
812	ENDDO
813	ENDDO
814	ENDDO
815	C
816	DO i = 1, klon
817	ztsol(i) = ftsol(i)
818	ENDDO
819	C
820	IF (if_ebil.ge.1) THEN
821	ztit='after dynamic'
822	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
823	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
824	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
825	C Comme les tendances de la physique sont ajoute dans la dynamique,
826	C on devrait avoir que la variation d'entalpie par la dynamique
827	C est egale a la variation de la physique au pas de temps precedent.
828	C Donc la somme de ces 2 variations devrait etre nulle.
829	call diagphy(cell_area,ztit,ip_ebil
830	e , zero_v, zero_v, zero_v, zero_v, zero_v
831	e , zero_v, zero_v, zero_v, ztsol
832	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
833	s , fs_bound, fq_bound )
834	END IF
835
836	c====================================================================
837	c XIOS outputs
838
839	#ifdef CPP_XIOS
840	! update XIOS time/calendar
841	call update_xios_timestep
842	#endif
843
844	c====================================================================
845	c Diagnostiquer la tendance dynamique
846	c
847	IF (ancien_ok) THEN
848	DO k = 1, klev
849	DO i = 1, klon
850	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
851	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
852	ENDDO
853	ENDDO
854
855	! ADAPTATION GCM POUR CP(T)
856	do i=1,klon
857	flux_dyn(i,1) = 0.0
858	do j=2,klev
859	flux_dyn(i,j) = flux_dyn(i,j-1)
860	. +cpnew(i,j-1)/RGd_t_dyn(i,j-1)(paprs(i,j-1)-paprs(i,j))
861	enddo
862	enddo
863
864	ELSE
865	DO k = 1, klev
866	DO i = 1, klon
867	d_u_dyn(i,k) = 0.0
868	d_t_dyn(i,k) = 0.0
869	ENDDO
870	ENDDO
871	ancien_ok = .TRUE.
872	ENDIF
873	c====================================================================
874
875	! Compute vertical velocity (Pa/s) from vertical mass flux
876	! Need to linearly interpolate mass flux to mid-layers
877	do k=1,klev-1
878	omega(1:klon,k) = 0.5RG(flxmw(1:klon,k)+flxmw(1:klon,k+1))
879	. / cell_area(1:klon)
880	enddo
881	omega(1:klon,klev) = 0.5RGflxmw(1:klon,klev) / cell_area(1:klon)
882
883	c======
884	c GEOP CORRECTION
885	c
886	c Ajouter le geopotentiel du sol:
887	c
888	DO k = 1, klev
889	DO i = 1, klon
890	zphi(i,k) = pphi(i,k) + pphis(i)
891	ENDDO
892	ENDDO
893
894	c............................
895	c CETTE CORRECTION VA DE PAIR AVEC DES MODIFS DE LEAPFROG(_p)
896	c ELLE MARCHE A 50 NIVEAUX (si mmean constante...)
897	c MAIS PAS A 78 NIVEAUX (quand mmean varie...)
898	c A ANALYSER PLUS EN DETAIL AVANT D'UTILISER
899	c............................
900	c zphi is recomputed (pphi is not ok if mean molecular mass varies)
901	c with dphi = RT/mmean d(ln p) [evaluated at interface]
902
903	c DO i = 1, klon
904	c zphi(i,1) = pphis(i) + Rt_seri(i,1)/mmean(i,1)1000.
905	c * *( log(paprs(i,1)) - log(pplay(i,1)) )
906	c DO k = 2, klev
907	c zphi(i,k) = zphi(i,k-1)
908	c * + R500.(t_seri(i,k)/mmean(i,k)+t_seri(i,k-1)/mmean(i,k-1))
909	c * * (log(pplay(i,k-1)) - log(pplay(i,k)))
910	c ENDDO
911	c ENDDO
912	c............................
913	c=====
914
915	c calcul de l'altitude aux niveaux intercouches
916	c ponderation des altitudes au niveau des couches en dp/p
917
918	DO i=1,klon
919	zzlay(i,1)=zphi(i,1)/RG ! [m]
920	zzlev(i,1)=pphis(i)/RG ! [m]
921	ENDDO
922	DO k=2,klev
923	DO i=1,klon
924	tlaymean=t_seri(i,k)
925	IF (t_seri(i,k).ne.t_seri(i,k-1))
926	& tlaymean=(t_seri(i,k)-t_seri(i,k-1))
927	& /log(t_seri(i,k)/t_seri(i,k-1))
928
929	zzlay(i,k)=zzlay(i,k-1)
930	& -(log(pplay(i,k)/pplay(i,k-1))rnew(i,k-1)tlaymean
931	& /(RG(RA/(RA+zzlay(i,k-1)))*2))
932	ENDDO
933	ENDDO
934	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
935	! Old : from geopotential. Problem when mu varies in the upper atm...
936	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
937	! DO k=1,klev
938	! DO i=1,klon
939	! zzlay(i,k)=zphi(i,k)/RG ! [m]
940	! ENDDO
941	! ENDDO
942	! DO i=1,klon
943	! zzlev(i,1)=pphis(i)/RG ! [m]
944	! ENDDO
945	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
946	DO k=2,klev
947	DO i=1,klon
948	z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
949	z2=(paprs(i,k) +pplay(i,k))/(paprs(i,k) -pplay(i,k))
950	zzlev(i,k)=(z1zzlay(i,k-1)+z2zzlay(i,k))/(z1+z2)
951	ENDDO
952	ENDDO
953	DO i=1,klon
954	zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
955	ENDDO
956
957	c====================================================================
958	c
959	c Verifier les temperatures
960	c
961	CALL hgardfou(t_seri,ftsol,'debutphy')
962
963	c====================================================================
964	c Orbite et eclairement
965	c=======================
966
967	c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0.
968	c donc pas de variations de Ls, ni de dist.
969	c La seule chose qui compte, c'est la rotation de la planete devant
970	c le Soleil...
971
972	zlongi = 0.0
973	dist = 0.72333 ! en UA
974
975	c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite
976	c a sa valeur, et prendre en compte leur evolution,
977	c il faudra refaire un orbite.F...
978	c CALL orbite(zlongi,dist)
979
980	IF (cycle_diurne) THEN
981	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
982	CALL zenang(zlongi,gmtime,zdtime,latitude_deg,longitude_deg,
983	& rmu0,fract)
984	ELSE
985	call mucorr(klon,zlongi,latitude_deg,rmu0,fract)
986	ENDIF
987
988	! print fraction of venus day
989
990	if (is_master) then
991	print*, 'gmtime = ', gmtime
992	end if
993
994	c======================================================================
995	c FIN INITIALISATIONS
996	c======================================================================
997	c======================================================================
998
999	c=======================================================
1000	! CONDUCTION and MOLECULAR VISCOSITY
1001	c=======================================================
1002
1003	d_t_conduc(:,:)=0.
1004	d_u_molvis(:,:)=0.
1005	d_v_molvis(:,:)=0.
1006
1007	IF (callthermos) THEN
1008
1009	tsurf(:)=t_seri(:,1)
1010	call conduction(klon, klev,pdtphys,
1011	$ pplay,paprs,t_seri,
1012	$ tsurf,zzlev,zzlay,d_t_conduc)
1013
1014	call molvis(klon, klev,pdtphys,
1015	$ pplay,paprs,t_seri,
1016	$ u,tsurf,zzlev,zzlay,d_u_molvis)
1017
1018	call molvis(klon, klev, pdtphys,
1019	$ pplay,paprs,t_seri,
1020	$ v,tsurf,zzlev,zzlay,d_v_molvis)
1021
1022	DO k=1,klev
1023	DO ig=1,klon
1024	t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K]
1025	u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! [m/s]
1026	v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! [m/s]
1027	ENDDO
1028	ENDDO
1029	ENDIF
1030
1031	c====================================================================
1032	c Compute mean mass, cp and R :
1033	c------------------------------------
1034
1035	if(callthermos) then
1036	call concentrations2(pplay,t_seri,tr_seri, nqmax)
1037	endif
1038
1039
1040	c=======================================================
1041	! CHEMISTRY AND MICROPHYSICS
1042	c=======================================================
1043
1044	if (iflag_trac.eq.1) then
1045
1046	if ( ok_aoa ) then
1047	call age_of_air(tr_seri(:,:,i_aoa),klon,klev,
1048	I itap,pdtphys,init_age,
1049	O age)
1050	end if
1051
1052	!====================================================================
1053	! Case 1: pseudo-chemistry with relaxation toward fixed profile
1054	!=========
1055	if (tr_scheme.eq.1) then
1056
1057	call phytrac_relax (debut,lafin,nqmax,
1058	I nlon,nlev,dtime,pplay,
1059	O tr_seri)
1060
1061	elseif (tr_scheme.eq.2) then
1062	!====================================================================
1063	! Case 2: surface emission
1064	! For the moment, inspired from Mars version
1065	! However, the variable 'source' could be used in physiq
1066	! so the call to phytrac_emiss could be to initialise it.
1067	!=========
1068
1069	call phytrac_emiss (debut,lafin,nqmax,
1070	I nlon,nlev,dtime,paprs,
1071	I latitude_deg,longitude_deg,
1072	O tr_seri)
1073
1074	else if (tr_scheme.eq.3) then
1075	!====================================================================
1076	! Case 3: Full chemistry and/or clouds.
1077	! routines are called every "chempas" physical timestep.
1078	!
1079	! if the physics is called 96000 times per venus day:
1080	!
1081	! nbapp_chem = 24000 => chempas = 4 => zctime = 420 s
1082	! nbapp_chem = 12000 => chempas = 8 => zctime = 840 s
1083	!=========
1084
1085
1086	chempas = nint(rday/pdtphys/nbapp_chem)
1087	zctime = dtime*real(chempas) ! chemical timestep
1088
1089	if (mod(itap,chempas) == 0) then ! <------- start of chemistry supercycling
1090
1091	! photochemistry and microphysics
1092
1093	call phytrac_chimie(debut,
1094	$ gmtime,
1095	$ nqmax,
1096	$ klon,
1097	$ latitude_deg,
1098	$ longitude_deg,
1099	$ zzlay,
1100	$ nlev,
1101	$ zctime,
1102	$ t_seri,
1103	$ pplay,
1104	$ tr_seri,
1105	$ d_tr_chem,
1106	$ iter,
1107	$ prod_tr,
1108	$ loss_tr,
1109	$ no_emission,
1110	$ o2_emission)
1111
1112	if (ok_sedim) then
1113	if (cl_scheme == 1) then
1114
1115	! sedimentation for simplified microphysics
1116
1117	#ifndef MESOSCALE
1118	call new_cloud_sedim(nb_mode,
1119	$ klon,
1120	$ nlev,
1121	$ zctime,
1122	$ pplay,
1123	$ paprs,
1124	$ t_seri,
1125	$ tr_seri,
1126	$ d_tr_chem,
1127	$ d_tr_sed(:,:,1:2),
1128	$ nqmax,
1129	$ Fsedim)
1130
1131	! test to avoid nans
1132
1133	do k = 1, klev
1134	do i = 1, klon
1135	if ((d_tr_sed(i,k,1) /= d_tr_sed(i,k,1)) .or.
1136	$ (d_tr_sed(i,k,2) /= d_tr_sed(i,k,2))) then
1137	print*,'sedim NaN PROBLEM'
1138	print*,'d_tr_sed Nan?',d_tr_sed(i,k,:)
1139	print*,'Temp',t_seri(i,k)
1140	print*,'lat-lon',i,'level',k,'zctime',zctime
1141	print*,'F_sed',Fsedim(i,k)
1142	d_tr_sed(i,k,:) = 0.
1143	end if
1144	end do
1145	end do
1146
1147	! tendency due to condensation and sedimentation
1148
1149	d_tr_sed(:,:,1:2) = d_tr_sed(:,:,1:2)/zctime
1150	Fsedim(:,1:klev) = Fsedim(:,1:klev)/zctime
1151	Fsedim(:,klev+1) = 0.
1152
1153	else if (cl_scheme == 2) then
1154
1155	! sedimentation for detailed microphysics
1156
1157	d_tr_sed(:,:,:) = 0.
1158
1159	do i = 1, klon
1160
1161	! mode 1
1162
1163	m0_mode1(:,1) = tr_seri(i,:,i_m0_mode1drop)
1164	m0_mode1(:,2) = tr_seri(i,:,i_m0_mode1ccn)
1165	m3_mode1(:,1) = tr_seri(i,:,i_m3_mode1sa)
1166	m3_mode1(:,2) = tr_seri(i,:,i_m3_mode1w)
1167	m3_mode1(:,3) = tr_seri(i,:,i_m3_mode1ccn)
1168
1169	call drop_sedimentation(zctime,klev,m0_mode1,m3_mode1,
1170	$ paprs(i,:),zzlev(i,:),
1171	$ zzlay(i,:),t_seri(i,:),1,
1172	$ d_ccn_sed(:,1),d_drop_sed,
1173	$ d_ccn_sed(:,2),d_liq_sed)
1174
1175	d_tr_sed(i,:,i_m0_mode1drop)= d_tr_sed(i,:,i_m0_mode1drop)
1176	$ + d_drop_sed
1177	d_tr_sed(i,:,i_m0_mode1ccn) = d_tr_sed(i,:,i_m0_mode1ccn)
1178	$ + d_ccn_sed(:,1)
1179	d_tr_sed(i,:,i_m3_mode1ccn) = d_tr_sed(i,:,i_m3_mode1ccn)
1180	$ + d_ccn_sed(:,2)
1181	d_tr_sed(i,:,i_m3_mode1sa) = d_tr_sed(i,:,i_m3_mode1sa)
1182	$ + d_liq_sed(:,1)
1183	d_tr_sed(i,:,i_m3_mode1w) = d_tr_sed(i,:,i_m3_mode1w)
1184	$ + d_liq_sed(:,2)
1185
1186	! mode 2
1187
1188	m0_mode2(:,1) = tr_seri(i,:,i_m0_mode2drop)
1189	m0_mode2(:,2) = tr_seri(i,:,i_m0_mode2ccn)
1190	m3_mode2(:,1) = tr_seri(i,:,i_m3_mode2sa)
1191	m3_mode2(:,2) = tr_seri(i,:,i_m3_mode2w)
1192	m3_mode2(:,3) = tr_seri(i,:,i_m3_mode2ccn)
1193
1194	call drop_sedimentation(zctime,klev,m0_mode2,m3_mode2,
1195	$ paprs(i,:),zzlev(i,:),
1196	$ zzlay(i,:),t_seri(i,:),2,
1197	$ d_ccn_sed(:,1),d_drop_sed,
1198	$ d_ccn_sed(:,2),d_liq_sed)
1199
1200	d_tr_sed(i,:,i_m0_mode2drop)= d_tr_sed(i,:,i_m0_mode2drop)
1201	$ + d_drop_sed
1202	d_tr_sed(i,:,i_m0_mode2ccn) = d_tr_sed(i,:,i_m0_mode2ccn)
1203	$ + d_ccn_sed(:,1)
1204	d_tr_sed(i,:,i_m3_mode2ccn) = d_tr_sed(i,:,i_m3_mode2ccn)
1205	$ + d_ccn_sed(:,2)
1206	d_tr_sed(i,:,i_m3_mode2sa) = d_tr_sed(i,:,i_m3_mode2sa)
1207	$ + d_liq_sed(:,1)
1208	d_tr_sed(i,:,i_m3_mode2w) = d_tr_sed(i,:,i_m3_mode2w)
1209	$ + d_liq_sed(:,2)
1210
1211	! aer
1212
1213	call aer_sedimentation(zctime,klev,
1214	$ tr_seri(i,:,i_m0_aer),
1215	$ tr_seri(i,:,i_m3_aer),
1216	$ paprs(i,:),zzlev(i,:),
1217	$ zzlay(i,:),t_seri(i,:),
1218	$ d_tr_sed(i,:,i_m0_aer),
1219	$ d_tr_sed(i,:,i_m3_aer),
1220	$ aer_flux)
1221
1222	end do
1223
1224	! tendency due to sedimentation
1225
1226	do iq = nqmax-nmicro+1,nqmax
1227	d_tr_sed(:,:,iq) = d_tr_sed(:,:,iq)/zctime
1228	end do
1229	#endif
1230	end if ! cl_scheme
1231
1232	! update gaseous tracers (chemistry)
1233
1234	do iq = 1, nqmax - nmicro
1235	tr_seri(:,:,iq) = max(tr_seri(:,:,iq)
1236	$ + d_tr_chem(:,:,iq)*zctime,1.e-30)
1237	end do
1238
1239	! update condensed tracers (condensation + sedimentation)
1240
1241	if (cl_scheme == 1) then
1242	tr_seri(:,:,i_h2so4liq) = max(tr_seri(:,:,i_h2so4liq)
1243	$ + d_tr_sed(:,:,1)*zctime, 1.e-30)
1244	tr_seri(:,:,i_h2oliq) = max(tr_seri(:,:,i_h2oliq)
1245	$ + d_tr_sed(:,:,2)*zctime, 1.e-30)
1246	else if (cl_scheme == 2) then
1247	do iq = nqmax-nmicro+1,nqmax
1248	tr_seri(:,:,iq) = max(tr_seri(:,:,iq)
1249	$ + d_tr_sed(:,:,iq)*zctime,1.e-30)
1250	end do
1251	end if ! cl_scheme
1252
1253	end if ! ok_sedim
1254	end if ! mod(itap,chempas) <------- end of chemistry supercycling
1255
1256	!=========
1257	! End Case 3: Full chemistry and/or clouds.
1258	!====================================================================
1259
1260	end if ! tr_scheme
1261	end if ! iflag_trac
1262
1263	c====================================================================
1264	c Appeler la diffusion verticale (programme de couche limite)
1265	c====================================================================
1266
1267	c-------------------------------
1268	c VENUS TEST: on ne tient pas compte des calculs de clmain mais on force
1269	c l'equilibre radiatif du sol
1270	if (.not. ok_clmain) then
1271	if (debut) then
1272	print*,"ATTENTION, CLMAIN SHUNTEE..."
1273	endif
1274
1275	DO i = 1, klon
1276	sens(i) = 0.0e0 ! flux de chaleur sensible au sol
1277	fder(i) = 0.0e0
1278	dlw(i) = 0.0e0
1279	ENDDO
1280
1281	c Incrementer la temperature du sol
1282	c
1283	DO i = 1, klon
1284	d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
1285	ftsol(i) = ftsol(i) + d_ts(i)
1286	do j=1,nsoilmx
1287	ftsoil(i,j)=ftsol(i)
1288	enddo
1289	ENDDO
1290
1291	c-------------------------------
1292	else
1293	c-------------------------------
1294
1295	fder = dlw
1296
1297	! ADAPTATION GCM POUR CP(T)
1298
1299	if (physideal) then
1300	CALL clmain_ideal(dtime,itap,
1301	e t_seri,u_seri,v_seri,
1302	e rmu0,
1303	e ftsol,
1304	$ ftsoil,
1305	$ paprs,pplay,ppk,radsol,falbe,
1306	e solsw, sollw, sollwdown, fder,
1307	e longitude_deg, latitude_deg, dx, dy,
1308	e debut, lafin,
1309	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
1310	s fluxt,fluxu,fluxv,cdragh,cdragm,
1311	s dsens,
1312	s ycoefh,yu1,yv1)
1313	else
1314	CALL clmain(dtime,itap,
1315	e t_seri,u_seri,v_seri,
1316	e rmu0,
1317	e ftsol,
1318	$ ftsoil,
1319	$ paprs,pplay,ppk,radsol,falbe,
1320	e solsw, sollw, sollwdown, fder,
1321	e longitude_deg, latitude_deg, dx, dy,
1322	& q2,
1323	e debut, lafin,
1324	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
1325	s fluxt,fluxu,fluxv,cdragh,cdragm,
1326	s dsens,
1327	s ycoefh,yu1,yv1)
1328	endif
1329
1330	CXXX Incrementation des flux
1331	DO i = 1, klon
1332	sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
1333	fder(i) = dlw(i) + dsens(i)
1334	ENDDO
1335	CXXX
1336	IF (.not. turb_resolved) then !True only for LES
1337	DO k = 1, klev
1338	DO i = 1, klon
1339	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1340	d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s
1341	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1342	d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s
1343	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1344	d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s
1345	ENDDO
1346	ENDDO
1347	ENDIF
1348	C TRACEURS
1349
1350	if (iflag_trac.eq.1) then
1351	DO k = 1, klev
1352	DO i = 1, klon
1353	delp(i,k) = paprs(i,k)-paprs(i,k+1)
1354	ENDDO
1355	ENDDO
1356
1357	DO iq=1, nqmax
1358
1359	CALL cltrac(dtime,ycoefh,t_seri,
1360	s tr_seri(:,:,iq),source(:,iq),
1361	e paprs, pplay,delp,
1362	s d_tr_vdf(:,:,iq))
1363
1364	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
1365	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
1366
1367	ENDDO !nqmax
1368
1369	endif
1370
1371	IF (if_ebil.ge.2) THEN
1372	ztit='after clmain'
1373	CALL diagetpq(cell_area,ztit,ip_ebil,2,1,dtime
1374	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1375	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1376	call diagphy(cell_area,ztit,ip_ebil
1377	e , zero_v, zero_v, zero_v, zero_v, sens
1378	e , zero_v, zero_v, zero_v, ztsol
1379	e , d_h_vcol, d_qt, d_ec
1380	s , fs_bound, fq_bound )
1381	END IF
1382	C
1383	c
1384	c Incrementer la temperature du sol
1385	c
1386	DO i = 1, klon
1387	ftsol(i) = ftsol(i) + d_ts(i)
1388	ENDDO
1389
1390	c Calculer la derive du flux infrarouge
1391	c
1392	DO i = 1, klon
1393	dlw(i) = - 4.0RSIGMAftsol(i)**3
1394	ENDDO
1395
1396	c-------------------------------
1397	endif ! fin du VENUS TEST
1398
1399	! tests: output tendencies
1400	! call writefield_phy('physiq_d_t_vdf',d_t_vdf,klev)
1401	! call writefield_phy('physiq_d_u_vdf',d_u_vdf,klev)
1402	! call writefield_phy('physiq_d_v_vdf',d_v_vdf,klev)
1403	! call writefield_phy('physiq_d_ts',d_ts,1)
1404
1405	c===================================================================
1406	c Convection seche
1407	c===================================================================
1408	c
1409	d_t_ajs(:,:)=0.
1410	d_u_ajs(:,:)=0.
1411	d_v_ajs(:,:)=0.
1412	d_tr_ajs(:,:,:)=0.
1413	c
1414	IF(prt_level>9)WRITE(lunout,*)
1415	. 'AVANT LA CONVECTION SECHE , iflag_ajs='
1416	s ,iflag_ajs
1417
1418	if(iflag_ajs.eq.0) then
1419	c Rien
1420	c ====
1421	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
1422
1423	else if(iflag_ajs.eq.1) then
1424
1425	c Ajustement sec
1426	c ==============
1427	IF(prt_level>9)WRITE(lunout,*)'ajsec'
1428
1429	! ADAPTATION GCM POUR CP(T)
1430	CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
1431	. tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
1432
1433	! ADAPTATION GCM POUR CP(T)
1434	do i=1,klon
1435	flux_ajs(i,1) = 0.0
1436	do j=2,klev
1437	flux_ajs(i,j) = flux_ajs(i,j-1)
1438	. + cpnew(i,j-1)/RG*d_t_ajs(i,j-1)/dtime
1439	. *(paprs(i,j-1)-paprs(i,j))
1440	enddo
1441	enddo
1442
1443	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
1444	d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s
1445	u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
1446	d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s
1447	v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
1448	d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s
1449
1450	if (iflag_trac.eq.1) then
1451	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
1452	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
1453	endif
1454	endif
1455
1456	! tests: output tendencies
1457	! call writefield_phy('physiq_d_t_ajs',d_t_ajs,klev)
1458	! call writefield_phy('physiq_d_u_ajs',d_u_ajs,klev)
1459	! call writefield_phy('physiq_d_v_ajs',d_v_ajs,klev)
1460	c
1461	IF (if_ebil.ge.2) THEN
1462	ztit='after dry_adjust'
1463	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
1464	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1465	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1466	call diagphy(cell_area,ztit,ip_ebil
1467	e , zero_v, zero_v, zero_v, zero_v, sens
1468	e , zero_v, zero_v, zero_v, ztsol
1469	e , d_h_vcol, d_qt, d_ec
1470	s , fs_bound, fq_bound )
1471	END IF
1472
1473	c====================================================================
1474	c RAYONNEMENT
1475	c====================================================================
1476	if (mod(itap,radpas) == 0) then
1477
1478	dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1479
1480	! update mmean
1481	if (ok_chem) then
1482	mmean(:,:) = 0.
1483	do iq = 1,nqmax - nmicro
1484	mmean(:,:) = mmean(:,:)+tr_seri(:,:,iq)/m_tr(iq)
1485	enddo
1486	mmean(:,:) = 1./mmean(:,:)
1487	rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K
1488	endif
1489
1490	cc---------------------------------------------
1491	if (callnlte .or. callthermos) then
1492	if (ok_chem) then
1493
1494	! nlte : use computed chemical species
1495
1496	co2vmr_gcm(:,:) = tr_seri(:,:,i_co2)*mmean(:,:)/m_tr(i_co2)
1497	covmr_gcm(:,:) = tr_seri(:,:,i_co)*mmean(:,:)/m_tr(i_co)
1498	ovmr_gcm(:,:) = tr_seri(:,:,i_o)*mmean(:,:)/m_tr(i_o)
1499	n2vmr_gcm(:,:) = tr_seri(:,:,i_n2)*mmean(:,:)/m_tr(i_n2)
1500
1501	else
1502
1503	! nlte : use hedin climatology
1504
1505	call compo_hedin83_mod(pplay,rmu0,
1506	$ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
1507	end if
1508	end if
1509
1510	c NLTE cooling from CO2 emission
1511	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1512
1513	IF(callnlte) THEN
1514	if(nltemodel.eq.0.or.nltemodel.eq.1) then
1515	c nltecool call not correct...
1516	c CALL nltecool(klon, klev, nqmax, pplay*9.869e-6, t_seri,
1517	c $ tr_seri, d_t_nlte)
1518	abort_message='nltemodel=0 or 1 should not be used...'
1519	call abort_physic(modname,abort_message,1)
1520	else if(nltemodel.eq.2) then
1521	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1522	!! HEDIN instead of compo for this calculation
1523	! call compo_hedin83_mod(pplay,rmu0,
1524	! $ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
1525	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1526	CALL nlte_tcool(klon,klev,pplay*9.869e-6,
1527	$ t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
1528	$ ovmr_gcm,d_t_nlte,ierr_nlte,varerr )
1529	if(ierr_nlte.gt.0) then
1530	write(,)
1531	$ 'WARNING: nlte_tcool output with error message',
1532	$ 'ierr_nlte=',ierr_nlte,'varerr=',varerr
1533	write(,)'I will continue anyway'
1534	endif
1535	endif
1536
1537	ELSE
1538
1539	d_t_nlte(:,:)=0.
1540
1541	ENDIF
1542
1543	c Find number of layers for LTE radiation calculations
1544
1545	IF(callnlte .or. callnirco2)
1546	$ CALL nlthermeq(klon, klev, paprs, pplay)
1547
1548	cc---------------------------------------------
1549	c LTE radiative transfert / solar / IR matrix
1550	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1551	if (physideal) then
1552	CALL radlwsw_newtoncool(presnivs,t_seri)
1553	else
1554	CALL radlwsw
1555	e (dist, rmu0, fract, zzlev,
1556	e paprs, pplay,ftsol, t_seri)
1557	endif
1558
1559	c ALBEDO VARIATIONS: test for Yeon Joo Lee
1560	c increment to increase it for 20 Vd => +80%
1561	c heat(:,:)=heat(:,:)(1.+0.80((rjourvrai-356)+gmtime)/20.)
1562	c or to decrease it for 20 Vd => 1/1.8
1563	c heat(:,:)=heat(:,:)/(1.+0.80*((rjourvrai-356)+gmtime)/20.)
1564
1565	c ------------ ALBEDO VARIATIONS: scenarios for VCD
1566	c shape of relative variation from Lee et al 2019 (Fig 13b)
1567	c between 57 km (4e4 Pa) and 72 km (2.5e3 Pa), peak at 67 km (6e3 Pa)
1568	c do j=1,klev
1569	c factAlb = 0.
1570	c if ((presnivs(j).gt.6.e3).and.(presnivs(j).lt.4.e4)) then
1571	c factAlb = (log(presnivs(j))-log(4.e4))/(log(6.e3)-log(4.e4))
1572	c elseif ((presnivs(j).lt.6.e3).and.(presnivs(j).gt.2.5e3)) then
1573	c factAlb = (log(presnivs(j))-log(2.5e3))/(log(6.e3)-log(2.5e3))
1574	c endif
1575	c Increase by 50% (Minimum albedo)
1576	c heat(:,j)=heat(:,j)(1+factAlb0.5)
1577	c Decrease by 30% (Maximum albedo)
1578	c heat(:,j)=heat(:,j)(1-factAlb0.3)
1579	c enddo
1580	c ------------ END ALBEDO VARIATIONS
1581
1582	cc---------------------------------------------
1583	c CO2 near infrared absorption
1584	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1585
1586	d_t_nirco2(:,:)=0.
1587	if (callnirco2) then
1588	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1589	!! HEDIN instead of compo for this calculation
1590	! call compo_hedin83_mod(pplay,rmu0,
1591	! $ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
1592	! tr_hedin(:,:,:)=tr_seri(:,:,:)
1593	! tr_hedin(:,:,i_co2)=co2vmr_gcm(:,:)/mmean(:,:)*m_tr(i_co2)
1594	! tr_hedin(:,:,i_co) = covmr_gcm(:,:)/mmean(:,:)*m_tr(i_co)
1595	! tr_hedin(:,:,i_o) = ovmr_gcm(:,:)/mmean(:,:)*m_tr(i_o)
1596	! tr_hedin(:,:,i_n2) = n2vmr_gcm(:,:)/mmean(:,:)*m_tr(i_n2)
1597	! call nirco2abs (klon, klev, pplay, dist, nqmax, tr_hedin,
1598	! . rmu0, fract, d_t_nirco2)
1599	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1600	call nirco2abs (klon, klev, pplay, dist, nqmax, tr_seri,
1601	. rmu0, fract, d_t_nirco2)
1602	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1603	endif
1604
1605
1606	cc---------------------------------------------
1607	c Net atmospheric radiative heating rate (K.s-1)
1608	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1609
1610	IF(callnlte.or.callnirco2) THEN
1611	CALL blendrad(klon, klev, pplay,heat,
1612	& cool, d_t_nirco2,d_t_nlte, dtsw, dtlw)
1613	ELSE
1614	dtsw(:,:)=heat(:,:)
1615	dtlw(:,:)=-1*cool(:,:)
1616	ENDIF
1617
1618	DO k=1,klev
1619	DO i=1,klon
1620	d_t_rad(i,k) = dtsw(i,k) + dtlw(i,k) ! K/s
1621	ENDDO
1622	ENDDO
1623
1624
1625	cc---------------------------------------------
1626	c EUV heating rate (K.s-1)
1627	c ~~~~~~~~~~~~~~~~~~~~~~~~
1628
1629	d_t_euv(:,:)=0.
1630
1631	IF (callthermos) THEN
1632
1633	call euvheat(klon, klev, nqmax, t_seri,paprs,pplay,zzlay,
1634	$ rmu0,dtimerad,gmtime,
1635	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1636	!! HEDIN instead of compo for this calculation
1637	!! cf nlte_tcool and nirco2abs above !!
1638	! $ tr_hedin, d_tr, d_t_euv )
1639	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1640	$ tr_seri, d_tr, d_t_euv )
1641	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1642
1643	DO k=1,klev
1644	DO ig=1,klon
1645	d_t_rad(ig,k)=d_t_rad(ig,k)+d_t_euv(ig,k)
1646	ENDDO
1647	ENDDO
1648
1649	ENDIF ! callthermos
1650
1651	ENDIF ! radpas
1652	c====================================================================
1653	c
1654	c Ajouter la tendance des rayonnements (tous les pas)
1655	c
1656	DO k = 1, klev
1657	DO i = 1, klon
1658	t_seri(i,k) = t_seri(i,k) + d_t_rad(i,k) * dtime
1659	ENDDO
1660	ENDDO
1661
1662	! increment physics counter
1663
1664	itap = itap + 1
1665	c====================================================================
1666
1667
1668	c==============================
1669	! -- MOLECULAR DIFFUSION ---
1670	c==============================
1671
1672	d_q_moldif(:,:,:)=0
1673
1674	IF (callthermos .and. ok_chem) THEN
1675
1676	call moldiff_red(klon, klev, nqmax,
1677	& pplay,paprs,t_seri, tr_seri, pdtphys,
1678	& d_t_euv,d_t_conduc,d_q_moldif)
1679
1680
1681	! --- update tendencies tracers ---
1682
1683	DO iq = 1, nqmax
1684	DO k=1,klev
1685	DO ig=1,klon
1686	tr_seri(ig,k,iq)= max(tr_seri(ig,k,iq)+
1687	& d_q_moldif(ig,k,iq)*dtime,1.e-30)
1688	ENDDO
1689	ENDDO
1690	ENDDO
1691
1692	ENDIF ! callthermos & ok_chem
1693
1694	c====================================================================
1695
1696	c==================================
1697	! -- ION AMBIPOLAR DIFFUSION ---
1698	c==================================
1699
1700	d_q_iondif(:,:,:)=0
1701
1702	IF (callthermos .and. ok_chem .and. ok_ionchem) THEN
1703	IF (ok_iondiff) THEN
1704
1705	call iondiff_red(klon,klev,nqmax,pplay,paprs,
1706	& t_seri,tr_seri,pphis,
1707	& gmtime,latitude_deg,longitude_deg,
1708	& pdtphys,d_t_euv,d_t_conduc,d_q_moldif,
1709	& d_q_iondif)
1710
1711	!write (,) 'TITI EST PASSE PAR LA'
1712	! --- update tendencies tracers ---
1713
1714	DO iq = 1, nqmax
1715	IF (type_tr(iq) .eq. 2) THEN
1716	DO k=1,klev
1717	DO ig=1,klon
1718	tr_seri(ig,k,iq)= max(tr_seri(ig,k,iq)+
1719	& d_q_iondif(ig,k,iq)*dtime,1.e-30)
1720	ENDDO
1721	ENDDO
1722	ENDIF
1723	ENDDO
1724	ENDIF ! ok_iondiff
1725	ENDIF ! callthermos & ok_chem & ok_ionchem
1726
1727	c====================================================================
1728	! tests: output tendencies
1729	! call writefield_phy('physiq_dtrad',dtrad,klev)
1730
1731	IF (if_ebil.ge.2) THEN
1732	ztit='after rad'
1733	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
1734	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1735	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1736	call diagphy(cell_area,ztit,ip_ebil
1737	e , topsw, toplw, solsw, sollw, zero_v
1738	e , zero_v, zero_v, zero_v, ztsol
1739	e , d_h_vcol, d_qt, d_ec
1740	s , fs_bound, fq_bound )
1741	END IF
1742	c
1743
1744	c====================================================================
1745	c Calcul des gravity waves FLOTT
1746	c====================================================================
1747	c
1748	c if (ok_orodr.or.ok_gw_nonoro) then
1749
1750	c CALCUL DE N2
1751	c UTILISE LA RELATION ENTRE N2 ET STABILITE
1752	c N2 = RG/T (dT/dz+RG/cp(T))
1753	c ET DONC EN N'UTILISE QUE LA TEMPERATURE, PAS teta.
1754
1755	do i=1,klon
1756	do k=2,klev
1757	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1758	enddo
1759	enddo
1760	do i=1,klon
1761	do k=2,klev
1762	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1763	zdtlev(i,k) = t_seri(i,k)-t_seri(i,k-1)
1764	zdzlev(i,k) = (zzlay(i,k)-zzlay(i,k-1))
1765	zn2(i,k) = RG/ztlev(i,k) * ( zdtlev(i,k)/zdzlev(i,k)
1766	. + RG/cpnew(i,k) )
1767	zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
1768	enddo
1769	zn2(i,1) = 1.e-12 ! securite
1770	enddo
1771
1772	c endif
1773
1774	c ----------------------------ORODRAG
1775	IF (ok_orodr) THEN
1776	c
1777	c selection des points pour lesquels le shema est actif:
1778	igwd=0
1779	DO i=1,klon
1780	itest(i)=0
1781	c IF ((zstd(i).gt.10.0)) THEN
1782	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1783	itest(i)=1
1784	igwd=igwd+1
1785	idx(igwd)=i
1786	ENDIF
1787	ENDDO
1788	c igwdim=MAX(1,igwd)
1789	c
1790	c A ADAPTER POUR VENUS!!! [ TN: c'est fait ! ]
1791	CALL drag_noro(klon,klev,dtime,paprs,pplay,pphi,zn2,
1792	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
1793	e igwd,idx,itest,
1794	e t_seri, u_seri, v_seri,
1795	s zulow, zvlow, zustrdr, zvstrdr,
1796	s d_t_oro, d_u_oro, d_v_oro,
1797	s zublstrdr,zvblstrdr,znlow,zeff,zbl,
1798	s ztau,tau0,knu2,kbreak)
1799
1800	c print*,"d_u_oro=",d_u_oro(klon/2,:)
1801	c ajout des tendances
1802	t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
1803	d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s
1804	u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
1805	d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s
1806	v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
1807	d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s
1808	c
1809	ELSE
1810	d_t_oro = 0.
1811	d_u_oro = 0.
1812	d_v_oro = 0.
1813	zustrdr = 0.
1814	zvstrdr = 0.
1815	zublstrdr = 0.
1816	zvblstrdr = 0.
1817	znlow = 0.
1818	zeff = 0.
1819	zbl = 0
1820	knu2 = 0
1821	kbreak = 0
1822	ztau = 0
1823	tau0 = 0.
1824	c
1825	ENDIF ! fin de test sur ok_orodr
1826	c
1827	! tests: output tendencies
1828	! call writefield_phy('physiq_d_t_oro',d_t_oro,klev)
1829	! call writefield_phy('physiq_d_u_oro',d_u_oro,klev)
1830	! call writefield_phy('physiq_d_v_oro',d_v_oro,klev)
1831
1832	c ----------------------------OROLIFT
1833	IF (ok_orolf) THEN
1834	print*,"ok_orolf NOT IMPLEMENTED !"
1835	stop
1836	c
1837	c selection des points pour lesquels le shema est actif:
1838	igwd=0
1839	DO i=1,klon
1840	itest(i)=0
1841	IF ((zpic(i)-zmea(i)).GT.100.) THEN
1842	itest(i)=1
1843	igwd=igwd+1
1844	idx(igwd)=i
1845	ENDIF
1846	ENDDO
1847	c igwdim=MAX(1,igwd)
1848	c
1849	c A ADAPTER POUR VENUS!!!
1850	c CALL lift_noro(klon,klev,dtime,paprs,pplay,
1851	c e latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval,
1852	c e igwd,idx,itest,
1853	c e t_seri, u_seri, v_seri,
1854	c s zulow, zvlow, zustrli, zvstrli,
1855	c s d_t_lif, d_u_lif, d_v_lif )
1856
1857	c
1858	c ajout des tendances
1859	t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
1860	d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s
1861	u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
1862	d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s
1863	v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
1864	d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s
1865	c
1866	ELSE
1867	d_t_lif = 0.
1868	d_u_lif = 0.
1869	d_v_lif = 0.
1870	zustrli = 0.
1871	zvstrli = 0.
1872	c
1873	ENDIF ! fin de test sur ok_orolf
1874
1875	c ---------------------------- NON-ORO GRAVITY WAVES
1876	IF(ok_gw_nonoro) then
1877
1878	! Obsolete
1879	! but used for VCD 1.1
1880	! call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
1881	! e t_seri, u_seri, v_seri, paprs(klon/2+1,:),
1882	! o zustrhi,zvstrhi,
1883	! o d_t_hin, d_u_hin, d_v_hin)
1884
1885	! New routine based on Generic
1886	! used after VCD 1.1, for VCD 2.0
1887	call nonoro_gwd_ran(klon,klev,dtime,pplay,zn2,presnivs,
1888	e t_seri, u_seri, v_seri,
1889	o zustrhi,zvstrhi,
1890	o d_t_hin, d_u_hin, d_v_hin)
1891
1892	c ajout des tendances
1893
1894	t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
1895	d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s
1896	u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
1897	d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s
1898	v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
1899	d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s
1900
1901	ELSE
1902	d_t_hin = 0.
1903	d_u_hin = 0.
1904	d_v_hin = 0.
1905	zustrhi = 0.
1906	zvstrhi = 0.
1907
1908	ENDIF ! fin de test sur ok_gw_nonoro
1909
1910	! tests: output tendencies
1911	! call writefield_phy('physiq_d_t_hin',d_t_hin,klev)
1912	! call writefield_phy('physiq_d_u_hin',d_u_hin,klev)
1913	! call writefield_phy('physiq_d_v_hin',d_v_hin,klev)
1914
1915	c====================================================================
1916	c Transport de ballons
1917	c====================================================================
1918	if (ballons.eq.1) then
1919	CALL ballon(30,pdtphys,rjourvrai,gmtime*RDAY,
1920	& latitude_deg,longitude_deg,
1921	c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM)
1922	C t,pplay,u,v,zphi) ! alt above planet average radius
1923	endif !ballons
1924
1925	c====================================================================
1926	c Bilan de mmt angulaire
1927	c====================================================================
1928	if (bilansmc.eq.1) then
1929	CMODDEB FLOTT
1930	C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
1931	C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
1932
1933	DO i = 1, klon
1934	zustrph(i)=0.
1935	zvstrph(i)=0.
1936	zustrcl(i)=0.
1937	zvstrcl(i)=0.
1938	ENDDO
1939	DO k = 1, klev
1940	DO i = 1, klon
1941	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
1942	c (paprs(i,k)-paprs(i,k+1))/rg
1943	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
1944	c (paprs(i,k)-paprs(i,k+1))/rg
1945	zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
1946	c (paprs(i,k)-paprs(i,k+1))/rg
1947	zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
1948	c (paprs(i,k)-paprs(i,k+1))/rg
1949	ENDDO
1950	ENDDO
1951
1952	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime*RDAY,
1953	C ra,rg,romega,
1954	C latitude_deg,longitude_deg,pphis,
1955	C zustrdr,zustrli,zustrcl,
1956	C zvstrdr,zvstrli,zvstrcl,
1957	C paprs,u,v)
1958
1959	CCMODFIN FLOTT
1960	endif !bilansmc
1961
1962	c====================================================================
1963	c====================================================================
1964	c Calculer le transport de l'eau et de l'energie (diagnostique)
1965	c
1966	c A REVOIR POUR VENUS...
1967	c
1968	c CALL transp (paprs,ftsol,
1969	c e t_seri, q_seri, u_seri, v_seri, zphi,
1970	c s ve, vq, ue, uq)
1971	c
1972	c====================================================================
1973	c+jld ec_conser
1974	DO k = 1, klev
1975	DO i = 1, klon
1976	d_t_ec(i,k)=0.5/cpnew(i,k)
1977	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
1978	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
1979	d_t_ec(i,k) = d_t_ec(i,k)/dtime
1980	END DO
1981	END DO
1982	do i=1,klon
1983	flux_ec(i,1) = 0.0
1984	do j=2,klev
1985	flux_ec(i,j) = flux_ec(i,j-1)
1986	. +cpnew(i,j-1)/RGd_t_ec(i,j-1)(paprs(i,j-1)-paprs(i,j))
1987	enddo
1988	enddo
1989
1990	c-jld ec_conser
1991	c====================================================================
1992	IF (if_ebil.ge.1) THEN
1993	ztit='after physic'
1994	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
1995	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1996	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1997	C Comme les tendances de la physique sont ajoute dans la dynamique,
1998	C on devrait avoir que la variation d'entalpie par la dynamique
1999	C est egale a la variation de la physique au pas de temps precedent.
2000	C Donc la somme de ces 2 variations devrait etre nulle.
2001	call diagphy(cell_area,ztit,ip_ebil
2002	e , topsw, toplw, solsw, sollw, sens
2003	e , zero_v, zero_v, zero_v, ztsol
2004	e , d_h_vcol, d_qt, d_ec
2005	s , fs_bound, fq_bound )
2006	C
2007	d_h_vcol_phy=d_h_vcol
2008	C
2009	END IF
2010	C
2011	c=======================================================================
2012	c SORTIES
2013	c=======================================================================
2014
2015	c Convertir les incrementations en tendances
2016	c
2017	DO k = 1, klev
2018	DO i = 1, klon
2019	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
2020	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
2021	d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
2022	ENDDO
2023	ENDDO
2024	c
2025	DO iq = 1, nqmax
2026	DO k = 1, klev
2027	DO i = 1, klon
2028	d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
2029	ENDDO
2030	ENDDO
2031	ENDDO
2032
2033	c mise à jour rho,mmean pour sorties
2034	if(callthermos) then
2035	call concentrations2(pplay,t_seri,tr_seri, nqmax)
2036	endif
2037
2038	c calcul vitesse verticale en m/s
2039	DO k = 1, klev
2040	DO i = 1, klon
2041	vertwind(i,k) = -omega(i,k)/(rho(i,k)*RG)
2042	END DO
2043	END DO
2044
2045	c------------------------
2046	c Calcul moment cinetique
2047	c------------------------
2048	c TEST VENUS...
2049	c mangtot = 0.0
2050	c DO k = 1, klev
2051	c DO i = 1, klon
2052	c mang(i,k) = RA*cos(latitude(i))
2053	c . (u_seri(i,k)+RAcos(latitude(i))*ROMEGA)
2054	c . cell_area(i)(paprs(i,k)-paprs(i,k+1))/RG
2055	c mangtot=mangtot+mang(i,k)
2056	c ENDDO
2057	c ENDDO
2058	c print*,"Moment cinetique total = ",mangtot
2059	c
2060	c------------------------
2061	c
2062	c Sauvegarder les valeurs de t et u a la fin de la physique:
2063	c
2064	DO k = 1, klev
2065	DO i = 1, klon
2066	u_ancien(i,k) = u_seri(i,k)
2067	t_ancien(i,k) = t_seri(i,k)
2068	ENDDO
2069	ENDDO
2070	c
2071	c=============================================================
2072	c Ecriture des sorties
2073	c=============================================================
2074	#ifndef MESOSCALE
2075	#ifdef CPP_IOIPSL
2076
2077	#ifdef histhf
2078	#include "write_histhf.h"
2079	#endif
2080
2081	#ifdef histday
2082	#include "write_histday.h"
2083	#endif
2084
2085	#ifdef histmth
2086	#include "write_histmth.h"
2087	#endif
2088
2089	#ifdef histins
2090	#include "write_histins.h"
2091	#endif
2092
2093	#endif
2094
2095	! XIOS outputs
2096	! This can be done ANYWHERE in the physics routines !
2097
2098	#ifdef CPP_XIOS
2099	! Send fields to XIOS: (NB these fields must also be defined as
2100	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
2101
2102	! 2D fields
2103
2104	CALL send_xios_field("phis",pphis)
2105	cell_area_out(:)=cell_area(:)
2106	if (is_north_pole_phy) cell_area_out(1)=cell_area(1)/nbp_lon
2107	if (is_south_pole_phy) cell_area_out(klon)=cell_area(klon)/nbp_lon
2108	CALL send_xios_field("aire",cell_area_out)
2109	CALL send_xios_field("tsol",ftsol)
2110	CALL send_xios_field("psol",paprs(:,1))
2111	CALL send_xios_field("cdragh",cdragh)
2112	CALL send_xios_field("cdragm",cdragm)
2113
2114	CALL send_xios_field("tops",topsw)
2115	CALL send_xios_field("topl",toplw)
2116	CALL send_xios_field("sols",solsw)
2117	CALL send_xios_field("soll",sollw)
2118
2119	! 3D fields
2120
2121	CALL send_xios_field("temp",t_seri)
2122	CALL send_xios_field("pres",pplay)
2123	CALL send_xios_field("geop",zphi)
2124	CALL send_xios_field("vitu",u_seri)
2125	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
2126	CALL send_xios_field("vitv",-1.*v_seri)
2127	CALL send_xios_field("vitw",omega)
2128	CALL send_xios_field("vitwz",vertwind)
2129	CALL send_xios_field("Kz",ycoefh)
2130	CALL send_xios_field("mmean",mmean)
2131	CALL send_xios_field("rho",rho)
2132	CALL send_xios_field("BV2",zn2)
2133
2134	CALL send_xios_field("dudyn",d_u_dyn)
2135	CALL send_xios_field("duvdf",d_u_vdf)
2136	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
2137	CALL send_xios_field("dvvdf",-1.*d_v_vdf)
2138	CALL send_xios_field("duajs",d_u_ajs)
2139	CALL send_xios_field("dugwo",d_u_oro)
2140	CALL send_xios_field("dugwno",d_u_hin)
2141	CALL send_xios_field("dvgwno",-1.*d_v_hin)
2142	CALL send_xios_field("dumolvis",d_u_molvis)
2143	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
2144	CALL send_xios_field("dvmolvis",-1.*d_v_molvis)
2145	CALL send_xios_field("dtdyn",d_t_dyn)
2146	CALL send_xios_field("dtphy",d_t)
2147	CALL send_xios_field("dtvdf",d_t_vdf)
2148	CALL send_xios_field("dtajs",d_t_ajs)
2149	CALL send_xios_field("dtswr",dtsw)
2150	CALL send_xios_field("dtswrNLTE",d_t_nirco2)
2151	CALL send_xios_field("dtswrLTE",heat)
2152	CALL send_xios_field("dtlwr",dtlw)
2153	CALL send_xios_field("dtlwrNLTE",d_t_nlte)
2154	CALL send_xios_field("dtlwrLTE",-1.*cool)
2155	CALL send_xios_field("dteuv",d_t_euv)
2156	CALL send_xios_field("dtcond",d_t_conduc)
2157	CALL send_xios_field("dtec",d_t_ec)
2158
2159	CALL send_xios_field("SWnet",swnet(:,1:klev))
2160	CALL send_xios_field("LWnet",lwnet(:,1:klev))
2161	CALL send_xios_field("fluxvdf",fluxt)
2162	CALL send_xios_field("fluxdyn",flux_dyn)
2163	CALL send_xios_field("fluxajs",flux_ajs)
2164	CALL send_xios_field("fluxec",flux_ec)
2165
2166	! when using tracers
2167
2168	if (iflag_trac == 1) then
2169
2170	if (ok_aoa) then
2171	call send_xios_field("age",age)
2172	call send_xios_field("aoa",tr_seri(:,:,i_aoa))
2173	endif
2174
2175	! production and destruction rate, cm-3.s-1
2176	! Beware of the context*.xml file !!
2177	if ((tr_scheme == 3) .and. (ok_chem)) THEN
2178	do iq = 1,nqmax - nmicro
2179	if ((iq.eq.i_o).or.(iq.eq.i_co)) THEN
2180	call send_xios_field("prod_"//tname(iq),
2181	$ prod_tr(:,:,iq))
2182	call send_xios_field("loss_"//tname(iq),
2183	$ loss_tr(:,:,iq))
2184	end if
2185	!if (iq.eq.i_o) then
2186	! call send_xios_field('prod_o', prod_tr(:,:,iq))
2187	! call send_xios_field('loss_o', loss_tr(:,:,iq))
2188	!end if
2189	!if (iq.eq.i_co) then
2190	! call send_xios_field('prod_co', prod_tr(:,:,iq))
2191	! call send_xios_field('loss_co', loss_tr(:,:,iq))
2192	!end if
2193	end do
2194	end if
2195
2196	! tracers in gas phase, volume mixing ratio
2197
2198	do iq = 1,nqmax - nmicro
2199	call send_xios_field(tname(iq),
2200	$ tr_seri(:,:,iq)*mmean(:,:)/m_tr(iq))
2201	end do
2202
2203	! tracers in gas phase, column densities
2204
2205	do iq = 1,nqmax - nmicro
2206	col_dens_tr(:,iq)=0.
2207	if (type_tr(iq).eq.1) THEN
2208	do k = 1, klev
2209	col_dens_tr(:,iq) = col_dens_tr(:,iq) +
2210	$ tr_seri(:,k,iq)* (paprs(:,k)-paprs(:,k+1)) / RG
2211	end do
2212	call send_xios_field("col_"//tname(iq),col_dens_tr(:,iq))
2213	end if
2214	end do
2215
2216	! tracers in liquid phase, volume mixing ratio
2217
2218	if ((tr_scheme == 3) .and. (cl_scheme == 1)) THEN
2219	call send_xios_field(tname(i_h2oliq),
2220	$ tr_seri(:,:,i_h2oliq)*mmean(:,:)/m_tr(i_h2oliq))
2221	call send_xios_field(tname(i_h2so4liq),
2222	$ tr_seri(:,:,i_h2so4liq)*mmean(:,:)/m_tr(i_h2so4liq))
2223	if (ok_sedim) then
2224	call send_xios_field("Fsedim",fsedim(:,1:klev))
2225	end if
2226	end if
2227
2228	! aeronomical emissions
2229
2230	call send_xios_field("no_emis",no_emission)
2231	call send_xios_field("o2_emis",o2_emission)
2232
2233	! chemical iterations
2234
2235	if (tr_scheme.eq.3) call send_xios_field("iter",real(iter))
2236
2237	end if
2238
2239	IF (callthermos .and. ok_chem) THEN
2240	CALL send_xios_field("d_qmoldifCO2",d_q_moldif(:,:,i_co2))
2241	CALL send_xios_field("d_qmoldifO3p",d_q_moldif(:,:,i_o))
2242	CALL send_xios_field("d_qmoldifN2",d_q_moldif(:,:,i_n2))
2243	ENDIF
2244
2245	!! DEBUG
2246	! if (is_master) print*,"itauphy=",itap
2247	! if (itap.eq.10) lafin=.true.
2248
2249	if (lafin.and.is_omp_master) then
2250	write(,) "physiq: call xios_context_finalize"
2251	call xios_context_finalize
2252	endif
2253
2254	#endif
2255	#else
2256	! Outputs MESOSCALE
2257	CALL allocate_comm_wrf(klon,klev)
2258	comm_HR_SW(1:klon,1:klev) = dtsw(1:klon,1:klev)
2259	comm_HR_LW(1:klon,1:klev) = dtlw(1:klon,1:klev)
2260	comm_DT_RAD(1:klon,1:klev) = d_t_rad(1:klon,1:klev)
2261	IF (turb_resolved) THEN
2262	open(17,file='hrdyn.txt',form='formatted',status='old')
2263	rewind(17)
2264	DO k=1,klev
2265	read(17,*) dt_dyn(k)
2266	ENDDO
2267	close(17)
2268
2269	do i=1,klon
2270	d_t(i,:)=d_t(i,:)+dt_dyn(:)
2271	comm_HR_DYN(i,:) = dt_dyn(:)
2272	enddo
2273	ELSE
2274	comm_HR_DYN(1:klon,1:klev) = d_t_dyn(1:klon,1:klev)
2275	comm_DT_VDF(1:klon,1:klev) = d_t_vdf(1:klon,1:klev)
2276	comm_DT_AJS(1:klon,1:klev) = d_t_ajs(1:klon,1:klev)
2277	ENDIF
2278	comm_DT(1:klon,1:klev)=d_t(1:klon,1:klev)
2279	#endif
2280
2281
2282	c====================================================================
2283	c Si c'est la fin, il faut conserver l'etat de redemarrage
2284	c====================================================================
2285	c
2286	IF (lafin) THEN
2287	itau_phy = itau_phy + itap
2288	CALL phyredem ("restartphy.nc")
2289
2290	c--------------FLOTT
2291	CMODEB LOTT
2292	C FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
2293	C DU BILAN DE MOMENT ANGULAIRE.
2294	if (bilansmc.eq.1) then
2295	write(,)'Fermeture de aaam_bud.out (FL Vous parle)'
2296	close(27)
2297	close(28)
2298	endif !bilansmc
2299	CMODFIN
2300	c-------------
2301	c--------------SLEBONNOIS
2302	C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
2303	C DES BALLONS
2304	if (ballons.eq.1) then
2305	write(,)'Fermeture des ballons*.out'
2306	close(30)
2307	close(31)
2308	close(32)
2309	close(33)
2310	close(34)
2311	endif !ballons
2312	c-------------
2313	ENDIF
2314
2315	END SUBROUTINE physiq
2316
2317	END MODULE physiq_mod
2318

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

Download in other formats:

Original Format