Context Navigation

physiq_mod.F @ 3061

Last change on this file since 3061 was 3035, checked in by amartinez, 23 months 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

Rev	Line
[3]	1	!
[1549]	2	! $Id: $
[3]	3	!
[1549]	4	MODULE physiq_mod
	5
	6	IMPLICIT NONE
	7
	8	CONTAINS
	9
[3]	10	SUBROUTINE physiq (nlon,nlev,nqmax,
	11	. debut,lafin,rjourvrai,gmtime,pdtphys,
[97]	12	. paprs,pplay,ppk,pphi,pphis,presnivs,
[3]	13	. u,v,t,qx,
[2125]	14	. flxmw,
[3]	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
[953]	40	c gmtime--input-R-fraction de la journee (0 a 1)
[3]	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)
[1301]	53	c flxmw---input-R-flux de masse vertical en kg/s
[3]	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======================================================================
[101]	61	USE ioipsl
[849]	62	! USE histcom ! not needed; histcom is included in ioipsl
[101]	63	use dimphy
[1545]	64	USE geometry_mod,only: longitude, latitude, ! in radians
	65	& longitude_deg,latitude_deg, ! in degrees
	66	& cell_area,dx,dy
[892]	67	USE phys_state_var_mod ! Variables sauvegardees de la physique
[1621]	68	USE cpdet_phy_mod, only: cpdet, t2tpot
[1310]	69	USE chemparam_mod
	70	USE conc
[2464]	71	USE param_v4_h
[1310]	72	USE compo_hedin83_mod2
[2486]	73	use radlwsw_newtoncool_mod, only: radlwsw_newtoncool
[1506]	74	! use ieee_arithmetic
[1524]	75	use time_phylmdz_mod, only: annee_ref, day_ref, itau_phy
[1530]	76	use mod_grid_phy_lmdz, only: nbp_lon
[1621]	77	use infotrac_phy, only: iflag_trac, tname, ttext
[1642]	78	use vertical_layers_mod, only: pseudoalt
[1723]	79	use turb_mod, only : sens, turb_resolved
[2580]	80	use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
[2259]	81	use sed_and_prod_mad, only: aer_sedimentation, drop_sedimentation
[2836]	82	use iono_h, only: temp_elect, temp_ion
[1642]	83	#ifdef CPP_XIOS
	84	use xios_output_mod, only: initialize_xios_output,
	85	& update_xios_timestep,
	86	& send_xios_field
[1682]	87	use wxios, only: wxios_context_init, xios_context_finalize
[1642]	88	#endif
[1723]	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,
[2193]	97	& is_south_pole_phy,
	98	& is_master
[1723]	99	#endif
[101]	100	IMPLICIT none
	101	c======================================================================
	102	c CLEFS CPP POUR LES IO
	103	c =====================
[1723]	104	#ifndef MESOSCALE
[101]	105	c#define histhf
	106	#define histday
	107	#define histmth
	108	#define histins
[1723]	109	#endif
[101]	110	c======================================================================
[3]	111	#include "dimsoil.h"
	112	#include "clesphys.h"
	113	#include "iniprint.h"
	114	#include "timerad.h"
[815]	115	#include "tabcontrol.h"
[1310]	116	#include "nirdata.h"
	117	#include "hedin.h"
[3]	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
[1301]	159	REAL flxmw(klon,klev)
[1665]	160
[3]	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
[2193]	179
	180	integer,save :: itap ! physics counter
[3]	181	REAL delp(klon,klev) ! epaisseur d'une couche
[2523]	182	REAL omega(klon,klev) ! vitesse verticale en Pa/s (+ downward)
	183	REAL vertwind(klon,klev) ! vitesse verticale en m/s (+ upward)
[3]	184
	185	INTEGER igwd,idx(klon),itest(klon)
[2193]	186
[3]	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)
[2047]	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)
[3]	197
	198	c Pour calcul GW drag oro et nonoro: CALCUL de N2:
[1665]	199	real zdtlev(klon,klev),zdzlev(klon,klev)
	200	real ztlev(klon,klev)
[3]	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
[105]	223	REAL,save,allocatable :: source(:,:)
[3]	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
[1723]	228	REAL dsens(klon) ! derivee chaleur sensible
[3]	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
[1305]	234	REAL Fsedim(klon,klev+1) ! Flux de sedimentation (kg.m-2)
[3]	235
	236	c======================================================================
	237	c
	238	c Declaration des procedures appelees
	239	c
[2836]	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
[3]	250	EXTERNAL printflag
	251	EXTERNAL zenang
	252	EXTERNAL diagetpq
	253	EXTERNAL conf_phys
	254	EXTERNAL diagphy
	255	EXTERNAL mucorr
[1310]	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
[1442]	267	EXTERNAL moldiff_red
	268
[3]	269	c
	270	c Variables locales
	271	c
	272	CXXX PB
[2836]	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
[3]	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
[2836]	281	REAL tmpout(klon,klev) ! [K/s]
[3]	282	c
	283	REAL dist, rmu0(klon), fract(klon)
	284	REAL zdtime, zlongi
	285	c
[1718]	286	INTEGER i, k, iq, ig, j, ll, ilon, ilat, ilev, isoil
[3]	287	c
	288	REAL zphi(klon,klev)
[1301]	289	REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
[2686]	290	REAL tlaymean ! valeur temporaire pour calculer zzlay
[1310]	291	real tsurf(klon)
[1301]	292
[1310]	293	c va avec nlte_tcool
	294	INTEGER ierr_nlte
	295	REAL varerr
	296
[2193]	297	! photochemistry
	298
	299	integer :: chempas
[3035]	300	real :: zctime
[2193]	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
[3]	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
[1310]	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]
[2836]	337	real d_u_molvis(klon,klev) ! [m/s] /s
	338	real d_v_molvis(klon,klev) ! [m/s] /s
[1310]	339
[1442]	340	c Tendencies due to molecular diffusion
	341	real d_q_moldif(klon,klev,nqmax)
	342
[3035]	343	c Tendencies due to ambipolar ion diffusion
	344	real d_q_iondif(klon,klev,nqmax)
	345
[3]	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
[1305]	365	REAL :: tr_seri(klon,klev,nqmax)
[2464]	366	REAL :: tr_hedin(klon,klev,nqmax)
[1305]	367	REAL :: d_tr(klon,klev,nqmax)
[2523]	368	c pour sorties
	369	REAL :: col_dens_tr(klon,nqmax)
[2580]	370	REAL,allocatable,save :: prod_tr(:,:,:)
	371	REAL,allocatable,save :: loss_tr(:,:,:)
[1305]	372
[3]	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
[97]	377	REAL zsto, zout, zsto1, zsto2, zero
[3]	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
[2580]	412	c ALBEDO VARIATIONS (VCD)
	413	REAL factAlb
[3]	414	c TEST VENUS...
	415	REAL mang(klon,klev) ! moment cinetique
	416	REAL mangtot ! moment cinetique total
	417
[1572]	418	c cell_area for outputs in hist*
	419	REAL cell_area_out(klon)
[1723]	420	#ifdef MESOSCALE
	421	REAL :: dt_dyn(klev)
	422	#endif
[3]	423	c Declaration des constantes et des fonctions thermodynamiques
	424	c
	425	#include "YOMCST.h"
	426
	427	c======================================================================
[2464]	428	c======================================================================
[3]	429	c INITIALISATIONS
[2464]	430	c======================================================================
[3]	431
	432	modname = 'physiq'
	433	ok_sync=.TRUE.
	434
[119]	435	bilansmc = 0
[3]	436	ballons = 0
[892]	437	! NE FONCTIONNENT PAS ENCORE EN PARALLELE !!!
[1723]	438	#ifndef MESOSCALE
[892]	439	if (is_parallel) then
	440	bilansmc = 0
	441	ballons = 0
	442	endif
[1723]	443	#endif
[3]	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
[2464]	455	c======================================================================
[3]	456	c PREMIER APPEL SEULEMENT
	457	c========================
	458	IF (debut) THEN
[105]	459	allocate(source(klon,nqmax))
[2580]	460	allocate(prod_tr(klon,klev,nqmax))
	461	allocate(loss_tr(klon,klev,nqmax))
[2795]	462	allocate(no_emission(klon,klev))
	463	allocate(o2_emission(klon,klev))
[3]	464
[1682]	465	#ifdef CPP_XIOS
	466	! Initialize XIOS context
	467	write(,) "physiq: call wxios_context_init"
	468	CALL wxios_context_init
	469	#endif
	470
[2135]	471	! The call to suphec is now done in iniphysiq_mod (interface)
	472	! CALL suphec ! initialiser constantes et parametres phys.
[3]	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
[2193]	482	call phys_state_var_init(nqmax)
[3]	483	c
[1310]	484	c Initialising Hedin model for upper atm
	485	c (to be revised when coupled to chemistry) :
	486	call conc_init
[1723]	487
[2193]	488	! initialise physics counter
	489
	490	itap = 0
	491
[1723]	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
[3]	516	c
	517	c Lecture startphy.nc :
	518	c
[892]	519	CALL phyetat0 ("startphy.nc")
[1718]	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
[1723]	528	#endif
[3]	529
[815]	530	c dtime est defini dans tabcontrol.h et lu dans startphy
[150]	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 '
[1719]	537	c call abort_physic(modname,abort_message,1)
[152]	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----------------
[150]	543	dtime=pdtphys
	544	ENDIF
	545
[2193]	546	radpas = NINT(RDAY/pdtphys/nbapp_rad)
[3]	547
[815]	548	CALL printflag( ok_journe,ok_instan )
[1642]	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
[3]	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
[1442]	599	source(:,:) = 0.0 ! pas de source, pour l'instant
[3]	600	c---------
[1310]	601
	602	c---------
	603	c INITIALIZE THERMOSPHERIC PARAMETERS
	604	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	605
	606	if (callthermos) then
[2464]	607	call fill_data_thermos
	608	call allocate_param_thermos(klev)
	609	call param_read_e107
[1310]	610	endif
	611
	612	c Initialisation (recomputed in concentration2)
	613	do ig=1,klon
	614	do j=1,klev
[2788]	615	rnew(ig,j)=RD
[1665]	616	cpnew(ig,j)=cpdet(t(ig,j))
[1442]	617	mmean(ig,j)=RMD
[1310]	618	akknew(ig,j)=1.e-4
[2788]	619	rho(ig,j)=pplay(ig,j)/(rnew(ig,j)*t(ig,j))
[1310]	620	enddo
	621
	622	enddo
	623
	624	IF(callthermos.or.callnlte.or.callnirco2) THEN
	625	call compo_hedin83_init2
	626	ENDIF
[1442]	627	if (callnlte.and.nltemodel.eq.2) call nlte_setup
	628	if (callnirco2.and.nircorr.eq.1) call nir_leedat
[1310]	629	c---------
	630
[3]	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'
[1719]	638	call abort_physic(modname,abort_message,1)
[3]	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'
[1719]	644	call abort_physic(modname,abort_message,1)
[3]	645	ENDIF
	646	c
[808]	647	IF (dtimeREAL(radpas).GT.(RDAY0.25).AND.cycle_diurne)
[3]	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'
[1719]	652	call abort_physic(modname,abort_message,1)
[3]	653	ENDIF
	654	c
	655	WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
	656	. iflag_ajs
	657	c
[1120]	658	ecrit_mth = NINT(RDAY/dtime*ecriphy) ! tous les ecritphy jours
[3]	659	IF (ok_mensuel) THEN
	660	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
	661	. ecrit_mth
	662	ENDIF
[97]	663
[3]	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
[97]	669
[1120]	670	ecrit_ins = NINT(RDAY/dtime*ecriphy) ! Fraction de jour reglable
[3]	671	IF (ok_instan) THEN
	672	WRITE(lunout,*)'La frequence de sortie instant. est de ',
	673	. ecrit_ins
	674	ENDIF
	675
[3035]	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
[3]	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
[2683]	719	! initialisation of microphysical and chemical parameters
[1661]	720
[2683]	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
[1661]	732
[2683]	733	! number of microphysical tracers
[2193]	734
	735	nmicro = 0
	736	if (ok_cloud .and. (cl_scheme == 1)) nmicro = 2
	737	if (ok_cloud .and. (cl_scheme == 2)) nmicro = 12
[1661]	738
[2683]	739	! initialise chemical parameters. includes the indexation of microphys tracers
[1305]	740
[2683]	741	if (ok_chem .or. cl_scheme == 2) then
	742	call chemparam_ini()
	743	end if
[3]	744
[2683]	745	! initialise cloud parameters (for cl_scheme = 1)
[1442]	746
[2683]	747	if (ok_cloud .and. (cl_scheme == 1)) then
	748	call cloud_ini(nlon,nlev)
	749	end if
	750
[2286]	751	! initialise mmean
	752
[2683]	753	if (callthermos) then
[2464]	754	call concentrations2(pplay,t,qx,nqmax)
[2683]	755	end if
[2286]	756
[2464]	757	c========================
[3]	758	ENDIF ! debut
[2464]	759	c========================
[2193]	760
[2464]	761	c======================================================================
[1642]	762	! ------------------------------------------------------
	763	! Initializations done at every physical timestep:
	764	! ------------------------------------------------------
[3]	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'
[1543]	809	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
[3]	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.
[1543]	816	call diagphy(cell_area,ztit,ip_ebil
[3]	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====================================================================
[1642]	824	c XIOS outputs
	825
	826	#ifdef CPP_XIOS
	827	! update XIOS time/calendar
	828	call update_xios_timestep
	829	#endif
	830
	831	c====================================================================
[3]	832	c Diagnostiquer la tendance dynamique
	833	c
	834	IF (ancien_ok) THEN
	835	DO k = 1, klev
[2836]	836	DO i = 1, klon
[3]	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
[2836]	839	ENDDO
[3]	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)
[2686]	847	. +cpnew(i,j-1)/RGd_t_dyn(i,j-1)(paprs(i,j-1)-paprs(i,j))
[3]	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====================================================================
[1301]	861
[2523]	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)
[1301]	869
[2047]	870	c======
	871	c GEOP CORRECTION
[3]	872	c
	873	c Ajouter le geopotentiel du sol:
	874	c
	875	DO k = 1, klev
[2836]	876	DO i = 1, klon
[3]	877	zphi(i,k) = pphi(i,k) + pphis(i)
[2836]	878	ENDDO
[3]	879	ENDDO
[1301]	880
[2047]	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
[2686]	902	c calcul de l'altitude aux niveaux intercouches
[1301]	903	c ponderation des altitudes au niveau des couches en dp/p
	904
	905	DO i=1,klon
[2686]	906	zzlay(i,1)=zphi(i,1)/RG ! [m]
[1310]	907	zzlev(i,1)=pphis(i)/RG ! [m]
[1301]	908	ENDDO
	909	DO k=2,klev
[2686]	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
[1301]	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
[3]	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
[2464]	952	c=======================
[3]	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
[2464]	960	dist = 0.72333 ! en UA
[3]	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
[808]	968	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
[1545]	969	CALL zenang(zlongi,gmtime,zdtime,latitude_deg,longitude_deg,
	970	& rmu0,fract)
[3]	971	ELSE
[1545]	972	call mucorr(klon,zlongi,latitude_deg,rmu0,fract)
[3]	973	ENDIF
	974
[2193]	975	! print fraction of venus day
[3]	976
[2193]	977	if (is_master) then
	978	print*, 'gmtime = ', gmtime
	979	end if
	980
[2464]	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]
[2836]	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]
[2464]	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
[3]	1031	if (iflag_trac.eq.1) then
[2193]	1032	!====================================================================
	1033	! Case 1: pseudo-chemistry with relaxation toward fixed profile
[2464]	1034	!=========
[1160]	1035	if (tr_scheme.eq.1) then
[1305]	1036
[1160]	1037	call phytrac_relax (debut,lafin,nqmax,
	1038	I nlon,nlev,dtime,pplay,
[3]	1039	O tr_seri)
[1160]	1040
	1041	elseif (tr_scheme.eq.2) then
[2193]	1042	!====================================================================
[1160]	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.
[2464]	1047	!=========
[1305]	1048
[2464]	1049	call phytrac_emiss (debut,lafin,nqmax,
[1160]	1050	I nlon,nlev,dtime,paprs,
[1545]	1051	I latitude_deg,longitude_deg,
[1160]	1052	O tr_seri)
[3]	1053
[2193]	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.
[2262]	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
[2464]	1063	!=========
[1305]	1064
[3035]	1065
[2193]	1066	chempas = nint(rday/pdtphys/nbapp_chem)
[2262]	1067	zctime = dtime*real(chempas) ! chemical timestep
[1305]	1068
[2193]	1069	if (mod(itap,chempas) == 0) then ! <------- start of chemistry supercycling
[1305]	1070
[2193]	1071	! photochemistry and microphysics
[1661]	1072
[2193]	1073	call phytrac_chimie(debut,
	1074	$ gmtime,
	1075	$ nqmax,
	1076	$ klon,
	1077	$ latitude_deg,
	1078	$ longitude_deg,
[2780]	1079	$ zzlay,
[2193]	1080	$ nlev,
	1081	$ zctime,
	1082	$ t_seri,
	1083	$ pplay,
	1084	$ tr_seri,
	1085	$ d_tr_chem,
[2580]	1086	$ iter,
	1087	$ prod_tr,
[2795]	1088	$ loss_tr,
	1089	$ no_emission,
	1090	$ o2_emission)
[2193]	1091
[1305]	1092	if (ok_sedim) then
[2193]	1093	if (cl_scheme == 1) then
[1305]	1094
[2193]	1095	! sedimentation for simplified microphysics
[1661]	1096
[1723]	1097	#ifndef MESOSCALE
[2193]	1098	call new_cloud_sedim(klon,
	1099	$ nlev,
	1100	$ zctime,
	1101	$ pplay,
	1102	$ paprs,
	1103	$ t_seri,
	1104	$ tr_seri,
[2200]	1105	$ d_tr_chem,
[2193]	1106	$ d_tr_sed(:,:,1:2),
	1107	$ nqmax,
	1108	$ Fsedim)
[1305]	1109
[2193]	1110	! test to avoid nans
[1442]	1111
[2193]	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
[1305]	1125
[2203]	1126	! tendency due to condensation and sedimentation
[1305]	1127
[2203]	1128	d_tr_sed(:,:,1:2) = d_tr_sed(:,:,1:2)/zctime
[2193]	1129	Fsedim(:,1:klev) = Fsedim(:,1:klev)/zctime
	1130	Fsedim(:,klev+1) = 0.
[1661]	1131
[2193]	1132	else if (cl_scheme == 2) then
[1661]	1133
[2193]	1134	! sedimentation for detailed microphysics
[1661]	1135
[2193]	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
[1661]	1154	d_tr_sed(i,:,i_m0_mode1drop)= d_tr_sed(i,:,i_m0_mode1drop)
[2193]	1155	$ + d_drop_sed
[1661]	1156	d_tr_sed(i,:,i_m0_mode1ccn) = d_tr_sed(i,:,i_m0_mode1ccn)
[2193]	1157	$ + d_ccn_sed(:,1)
[1661]	1158	d_tr_sed(i,:,i_m3_mode1ccn) = d_tr_sed(i,:,i_m3_mode1ccn)
[2193]	1159	$ + d_ccn_sed(:,2)
[1661]	1160	d_tr_sed(i,:,i_m3_mode1sa) = d_tr_sed(i,:,i_m3_mode1sa)
[2193]	1161	$ + d_liq_sed(:,1)
[1661]	1162	d_tr_sed(i,:,i_m3_mode1w) = d_tr_sed(i,:,i_m3_mode1w)
[2193]	1163	$ + d_liq_sed(:,2)
[1661]	1164
[2193]	1165	! mode 2
[1661]	1166
[2193]	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
[1661]	1179	d_tr_sed(i,:,i_m0_mode2drop)= d_tr_sed(i,:,i_m0_mode2drop)
[2193]	1180	$ + d_drop_sed
[1661]	1181	d_tr_sed(i,:,i_m0_mode2ccn) = d_tr_sed(i,:,i_m0_mode2ccn)
[2193]	1182	$ + d_ccn_sed(:,1)
[1661]	1183	d_tr_sed(i,:,i_m3_mode2ccn) = d_tr_sed(i,:,i_m3_mode2ccn)
[2193]	1184	$ + d_ccn_sed(:,2)
[1661]	1185	d_tr_sed(i,:,i_m3_mode2sa) = d_tr_sed(i,:,i_m3_mode2sa)
[2193]	1186	$ + d_liq_sed(:,1)
[1661]	1187	d_tr_sed(i,:,i_m3_mode2w) = d_tr_sed(i,:,i_m3_mode2w)
[2193]	1188	$ + d_liq_sed(:,2)
[1661]	1189
[2193]	1190	! aer
[1661]	1191
[2193]	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
[1661]	1202
[2193]	1203	! tendency due to sedimentation
	1204
[2194]	1205	do iq = nqmax-nmicro+1,nqmax
[2193]	1206	d_tr_sed(:,:,iq) = d_tr_sed(:,:,iq)/zctime
	1207	end do
[1723]	1208	#endif
[2193]	1209	end if ! cl_scheme
[1661]	1210
[2203]	1211	! update gaseous tracers (chemistry)
[1661]	1212
[2203]	1213	do iq = 1, nqmax - nmicro
[2464]	1214	tr_seri(:,:,iq) = max(tr_seri(:,:,iq)
	1215	$ + d_tr_chem(:,:,iq)*zctime,1.e-30)
[2203]	1216	end do
[1305]	1217
[2203]	1218	! update condensed tracers (condensation + sedimentation)
[1305]	1219
[2194]	1220	if (cl_scheme == 1) then
[2203]	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)
[2194]	1225	else if (cl_scheme == 2) then
[2203]	1226	do iq = nqmax-nmicro+1,nqmax
[2464]	1227	tr_seri(:,:,iq) = max(tr_seri(:,:,iq)
	1228	$ + d_tr_sed(:,:,iq)*zctime,1.e-30)
[2203]	1229	end do
[2194]	1230	end if ! cl_scheme
[2203]	1231
[2194]	1232	end if ! ok_sedim
[2203]	1233	end if ! mod(itap,chempas) <------- end of chemistry supercycling
	1234
[2464]	1235	!=========
[2194]	1236	! End Case 3: Full chemistry and/or clouds.
	1237	!====================================================================
[2193]	1238
	1239	end if ! tr_scheme
	1240	end if ! iflag_trac
	1241
[3]	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
[2135]	1249	if (.not. ok_clmain) then
[3]	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)
[808]	1277
[2135]	1278	if (physideal) then
	1279	CALL clmain_ideal(dtime,itap,
[3]	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,
[1545]	1286	e longitude_deg, latitude_deg, dx, dy,
[3]	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)
[2135]	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,
[2535]	1300	e longitude_deg, latitude_deg, dx, dy,
	1301	& q2,
[2135]	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
[3]	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
[1723]	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
[3]	1327	C TRACEURS
	1328
	1329	if (iflag_trac.eq.1) then
[101]	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
[1442]	1335
[3]	1336	DO iq=1, nqmax
[1442]	1337
[3]	1338	CALL cltrac(dtime,ycoefh,t_seri,
[1442]	1339	s tr_seri(:,:,iq),source(:,iq),
[3]	1340	e paprs, pplay,delp,
[1442]	1341	s d_tr_vdf(:,:,iq))
	1342
[3]	1343	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
	1344	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
[1442]	1345
	1346	ENDDO !nqmax
[3]	1347
[1442]	1348	endif
	1349
[3]	1350	IF (if_ebil.ge.2) THEN
	1351	ztit='after clmain'
[1543]	1352	CALL diagetpq(cell_area,ztit,ip_ebil,2,1,dtime
[3]	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)
[1543]	1355	call diagphy(cell_area,ztit,ip_ebil
[3]	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
[973]	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
[3]	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)
[2686]	1417	. + cpnew(i,j-1)/RG*d_t_ajs(i,j-1)/dtime
[3]	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
[1301]	1428
	1429	if (iflag_trac.eq.1) then
[3]	1430	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
	1431	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
[1301]	1432	endif
[3]	1433	endif
[973]	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)
[3]	1439	c
	1440	IF (if_ebil.ge.2) THEN
	1441	ztit='after dry_adjust'
[1543]	1442	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
[3]	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)
[1543]	1445	call diagphy(cell_area,ztit,ip_ebil
[3]	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
[2193]	1452	c====================================================================
	1453	c RAYONNEMENT
	1454	c====================================================================
	1455	if (mod(itap,radpas) == 0) then
[1310]	1456
[2193]	1457	dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s)
[1310]	1458
[2464]	1459	! update mmean
	1460	if (ok_chem) then
	1461	mmean(:,:) = 0.
	1462	do iq = 1,nqmax - nmicro
[2622]	1463	mmean(:,:) = mmean(:,:)+tr_seri(:,:,iq)/m_tr(iq)
[2464]	1464	enddo
[2622]	1465	mmean(:,:) = 1./mmean(:,:)
	1466	rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K
[2464]	1467	endif
	1468
	1469	cc---------------------------------------------
[2193]	1470	if (callnlte .or. callthermos) then
	1471	if (ok_chem) then
[1442]	1472
[2193]	1473	! nlte : use computed chemical species
[1442]	1474
[2193]	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)
[1442]	1479
[2193]	1480	else
[1442]	1481
[2193]	1482	! nlte : use hedin climatology
[1442]	1483
[2193]	1484	call compo_hedin83_mod(pplay,rmu0,
	1485	$ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
	1486	end if
	1487	end if
[1442]	1488
	1489	c NLTE cooling from CO2 emission
	1490	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1491
[1310]	1492	IF(callnlte) THEN
	1493	if(nltemodel.eq.0.or.nltemodel.eq.1) then
[2464]	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)
[1310]	1499	else if(nltemodel.eq.2) then
[2464]	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	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1442]	1505	CALL nlte_tcool(klon,klev,pplay*9.869e-6,
[1310]	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
[1543]	1520	ENDIF
[1310]	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
[2464]	1527	cc---------------------------------------------
[1442]	1528	c LTE radiative transfert / solar / IR matrix
	1529	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2135]	1530	if (physideal) then
[2486]	1531	CALL radlwsw_newtoncool(presnivs,t_seri)
[2135]	1532	else
	1533	CALL radlwsw
	1534	e (dist, rmu0, fract, zzlev,
	1535	e paprs, pplay,ftsol, t_seri)
	1536	endif
[808]	1537
[2580]	1538	c ALBEDO VARIATIONS: test for Yeon Joo Lee
[2135]	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.)
[1301]	1543
[2580]	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
[2464]	1561	cc---------------------------------------------
[1442]	1562	c CO2 near infrared absorption
	1563	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1310]	1564
	1565	d_t_nirco2(:,:)=0.
	1566	if (callnirco2) then
[2464]	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	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1442]	1579	call nirco2abs (klon, klev, pplay, dist, nqmax, tr_seri,
	1580	. rmu0, fract, d_t_nirco2)
[2464]	1581	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1310]	1582	endif
	1583
	1584
[2464]	1585	cc---------------------------------------------
[1310]	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
[1543]	1593	dtsw(:,:)=heat(:,:)
	1594	dtlw(:,:)=-1*cool(:,:)
	1595	ENDIF
[1310]	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
[1442]	1612	call euvheat(klon, klev, nqmax, t_seri,paprs,pplay,zzlay,
[2464]	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	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1442]	1619	$ tr_seri, d_tr, d_t_euv )
[2464]	1620	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1543]	1621
[1310]	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====================================================================
[3]	1632	c
	1633	c Ajouter la tendance des rayonnements (tous les pas)
	1634	c
	1635	DO k = 1, klev
	1636	DO i = 1, klon
[1310]	1637	t_seri(i,k) = t_seri(i,k) + d_t_rad(i,k) * dtime
[3]	1638	ENDDO
	1639	ENDDO
[973]	1640
[2193]	1641	! increment physics counter
	1642
	1643	itap = itap + 1
[2464]	1644	c====================================================================
[2193]	1645
[1310]	1646
[2464]	1647	c==============================
[1442]	1648	! -- MOLECULAR DIFFUSION ---
[2464]	1649	c==============================
[1442]	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,
[2464]	1657	& d_t_euv,d_t_conduc,d_q_moldif)
[1442]	1658
	1659
	1660	! --- update tendencies tracers ---
	1661
	1662	DO iq = 1, nqmax
	1663	DO k=1,klev
	1664	DO ig=1,klon
[2464]	1665	tr_seri(ig,k,iq)= max(tr_seri(ig,k,iq)+
	1666	& d_q_moldif(ig,k,iq)*dtime,1.e-30)
[1442]	1667	ENDDO
	1668	ENDDO
	1669	ENDDO
	1670
	1671	ENDIF ! callthermos & ok_chem
	1672
[1310]	1673	c====================================================================
[3035]	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====================================================================
[973]	1707	! tests: output tendencies
	1708	! call writefield_phy('physiq_dtrad',dtrad,klev)
[3]	1709
	1710	IF (if_ebil.ge.2) THEN
	1711	ztit='after rad'
[1543]	1712	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
[3]	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)
[1543]	1715	call diagphy(cell_area,ztit,ip_ebil
[3]	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
[2047]	1727	c if (ok_orodr.or.ok_gw_nonoro) then
[1665]	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
[3]	1734	do i=1,klon
	1735	do k=2,klev
[1543]	1736	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
	1737	enddo
[3]	1738	enddo
	1739	do i=1,klon
	1740	do k=2,klev
[1665]	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)
[2686]	1743	zdzlev(i,k) = (zzlay(i,k)-zzlay(i,k-1))
[1665]	1744	zn2(i,k) = RG/ztlev(i,k) * ( zdtlev(i,k)/zdzlev(i,k)
[2686]	1745	. + RG/cpnew(i,k) )
[3]	1746	zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
[1543]	1747	enddo
[808]	1748	zn2(i,1) = 1.e-12 ! securite
[3]	1749	enddo
	1750
[2047]	1751	c endif
[3]	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
[2047]	1769	c A ADAPTER POUR VENUS!!! [ TN: c'est fait ! ]
[1306]	1770	CALL drag_noro(klon,klev,dtime,paprs,pplay,pphi,zn2,
[3]	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,
[2047]	1775	s d_t_oro, d_u_oro, d_v_oro,
	1776	s zublstrdr,zvblstrdr,znlow,zeff,zbl,
	1777	s ztau,tau0,knu2,kbreak)
[3]	1778
[808]	1779	c print*,"d_u_oro=",d_u_oro(klon/2,:)
[3]	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.
[1543]	1792	zustrdr = 0.
	1793	zvstrdr = 0.
[2047]	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.
[3]	1803	c
	1804	ENDIF ! fin de test sur ok_orodr
	1805	c
[973]	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
[3]	1811	c ----------------------------OROLIFT
	1812	IF (ok_orolf) THEN
[808]	1813	print*,"ok_orolf NOT IMPLEMENTED !"
	1814	stop
[3]	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,
[1545]	1830	c e latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval,
[3]	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
[2580]	1857	! Obsolete
	1858	! but used for VCD 1.1
[2622]	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,
[3]	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
[973]	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
[3]	1894	c====================================================================
	1895	c Transport de ballons
	1896	c====================================================================
	1897	if (ballons.eq.1) then
[1545]	1898	CALL ballon(30,pdtphys,rjourvrai,gmtime*RDAY,
	1899	& latitude_deg,longitude_deg,
[1543]	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
[3]	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
[953]	1931	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime*RDAY,
[3]	1932	C ra,rg,romega,
[1545]	1933	C latitude_deg,longitude_deg,pphis,
[3]	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
[2686]	1955	d_t_ec(i,k)=0.5/cpnew(i,k)
[3]	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)
[2686]	1965	. +cpnew(i,j-1)/RGd_t_ec(i,j-1)(paprs(i,j-1)-paprs(i,j))
[3]	1966	enddo
	1967	enddo
	1968
	1969	c-jld ec_conser
	1970	c====================================================================
	1971	IF (if_ebil.ge.1) THEN
	1972	ztit='after physic'
[1543]	1973	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
[3]	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.
[1543]	1980	call diagphy(cell_area,ztit,ip_ebil
[3]	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
[2464]	2012	c mise à jour rho,mmean pour sorties
	2013	if(callthermos) then
	2014	call concentrations2(pplay,t_seri,tr_seri, nqmax)
	2015	endif
	2016
[2523]	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
[3]	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
[1545]	2031	c mang(i,k) = RA*cos(latitude(i))
	2032	c . (u_seri(i,k)+RAcos(latitude(i))*ROMEGA)
[1543]	2033	c . cell_area(i)(paprs(i,k)-paprs(i,k+1))/RG
[3]	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=============================================================
[1723]	2053	#ifndef MESOSCALE
[3]	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
[1642]	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
[1665]	2081	! 2D fields
	2082
	2083	CALL send_xios_field("phis",pphis)
[1642]	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))
[1665]	2090	CALL send_xios_field("cdragh",cdragh)
	2091	CALL send_xios_field("cdragm",cdragm)
[1642]	2092
[1665]	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
[1642]	2100	CALL send_xios_field("temp",t_seri)
[1665]	2101	CALL send_xios_field("pres",pplay)
	2102	CALL send_xios_field("geop",zphi)
[1642]	2103	CALL send_xios_field("vitu",u_seri)
	2104	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
	2105	CALL send_xios_field("vitv",-1.*v_seri)
[1665]	2106	CALL send_xios_field("vitw",omega)
[2523]	2107	CALL send_xios_field("vitwz",vertwind)
[1665]	2108	CALL send_xios_field("Kz",ycoefh)
	2109	CALL send_xios_field("mmean",mmean)
	2110	CALL send_xios_field("rho",rho)
[2047]	2111	CALL send_xios_field("BV2",zn2)
[1665]	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)
[2836]	2120	CALL send_xios_field("dvgwno",-1.*d_v_hin)
[1665]	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
[1687]	2138	CALL send_xios_field("SWnet",swnet(:,1:klev))
	2139	CALL send_xios_field("LWnet",lwnet(:,1:klev))
[1665]	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
[2193]	2145	! when using tracers
[1687]	2146
[2193]	2147	if (iflag_trac == 1) then
[1665]	2148
[2580]	2149	! production and destruction rate, cm-3.s-1
	2150	! Beware of the context*.xml file !!
[2836]	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
[2580]	2169
[2193]	2170	! tracers in gas phase, volume mixing ratio
	2171
	2172	do iq = 1,nqmax - nmicro
	2173	call send_xios_field(tname(iq),
[2464]	2174	$ tr_seri(:,:,iq)*mmean(:,:)/m_tr(iq))
[2193]	2175	end do
	2176
[2523]	2177	! tracers in gas phase, column densities
[2795]	2178
[2523]	2179	do iq = 1,nqmax - nmicro
	2180	col_dens_tr(:,iq)=0.
[2836]	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
[2523]	2188	end do
	2189
[2193]	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),
[2464]	2194	$ tr_seri(:,:,i_h2oliq)*mmean(:,:)/m_tr(i_h2oliq))
[2193]	2195	call send_xios_field(tname(i_h2so4liq),
[2464]	2196	$ tr_seri(:,:,i_h2so4liq)*mmean(:,:)/m_tr(i_h2so4liq))
[2193]	2197	if (ok_sedim) then
	2198	call send_xios_field("Fsedim",fsedim(:,1:klev))
	2199	end if
	2200	end if
	2201
[2795]	2202	! aeronomical emissions
	2203
[2836]	2204	call send_xios_field("no_emis",no_emission)
	2205	call send_xios_field("o2_emis",o2_emission)
[2795]	2206
[2193]	2207	! chemical iterations
	2208
[2464]	2209	if (tr_scheme.eq.3) call send_xios_field("iter",real(iter))
[2193]	2210
	2211	end if
	2212
[1665]	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
[1682]	2218
[2464]	2219	!! DEBUG
	2220	! if (is_master) print*,"itauphy=",itap
	2221	! if (itap.eq.10) lafin=.true.
	2222
[1682]	2223	if (lafin.and.is_omp_master) then
	2224	write(,) "physiq: call xios_context_finalize"
	2225	call xios_context_finalize
	2226	endif
	2227
[1642]	2228	#endif
[1723]	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)
[1642]	2242
[1723]	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
[3]	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
[892]	2262	CALL phyredem ("restartphy.nc")
[3]	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
[1549]	2289	END SUBROUTINE physiq
[3]	2290
[1549]	2291	END MODULE physiq_mod
[3]	2292

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

Download in other formats:

Original Format