Context Navigation

physiq_mod.F @ 3221

Last change on this file since 3221 was 3035, checked in by amartinez, 2 years ago

======= VENUS PCM ===============

COMMIT BY ANTOINE MARTINEZ

Implementation of vertical ambipolar diffusion in physics

=================================

NEW KEYWORD OF PHYSIQ.DEF

=================================

NEW VERSION OF "physiq.def"

deftank/physiq-96x96x90-chemistry-IONOSPHERE-IONDIFF.def

ok_iondiff: keyword supposed to activate ion ambipolar diffusion

nbapp_chem: replaces nbapp_chim, in order to characterize the Number of calls to the chemistry routines (per Venusian day)

================

phyvenus

================

Iondiff_red.F

Calculation of the Ion Ambipolar Diffusion for 13 ions on 14:

O2+, O+, C+, N+, CO2+,
NO+, CO+, H+, H2O+, H3O+,

HCO+, N2+, OH+

The ion temperature is fixed as the half of the electron temperature
calculated in ion_h.F for the stability of the program and because the
ion temperature is lower than electron temperature.

plasmaphys_venus_mod.F90

parameters of the ambipolar diffusion scheme:

parameter (Pdiffion = 7.0D-4) ! pressure in Pa below which ion diffusion is computed
parameter (naltvert = 168) ! number of level on the altimetric subgrid
parameter (nsubvert = 24000) ! ptimephysiq/nsubvert - minimum sub-timestep allowed
parameter ( nsubmin = 40) ! ptimephysiq/nsubmin - maximum sub-timestep allowed

physic_mod.F

nbapp_chem is not fixed anymore here but deftank/in physiq.def
Ambipolar diffusion activated if (ok_iondiff) is true

conf_phys.f90

add ok_iondiff as parameters to read from physiq.def file set to .false. by default
add nbapp_chem as parameters to characterize the Number of calls to the chemistry routines (per Venusian day) instead of be fixed in physic_mod.F

to read from physiq.def file set to 24000 by default

cleshphys.h

added ok_iondiff & nbapp_chem in COMMON/clesphys_l/
removed nbapp_chim

phytrac_chemistry.F

Added security in the calculation of the sza_local in order to avoid the rare case where the |range| is above 1

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

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

Download in other formats:

Original Format