Context Navigation

physiq_mod.F @ 2203

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

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

Download in other formats:

Original Format