Context Navigation

physiq_mod.F @ 2665

Last change on this file since 2665 was 2622, checked in by slebonnois, 3 years ago
SL: VENUS update (i) bug correction (2 bugs, phytrac and physiq), affected meam molec mass computations... (ii) updates for VCD 2.0 (iii) aeropacity: for latitudinal variations of the cloud distribution
File size: 70.9 KB

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

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

Download in other formats:

Original Format