Context Navigation

physiq_mod.F @ 2200

Last change on this file since 2200 was 2200, checked in by flefevre, 5 years ago
correction de bugs dans le supercycling de la chimie
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	REAL :: d_tr_ssed(klon)
303	c
304	c Variables du changement
305	c
306	c ajs: ajustement sec
307	c vdf: couche limite (Vertical DiFfusion)
308	REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax)
309	REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
310	c
311	REAL d_ts(klon)
312	c
313	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
314	REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax)
315	c
316	CMOD LOTT: Tendances Orography Sous-maille
317	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
318	REAL d_t_oro(klon,klev)
319	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
320	REAL d_t_lif(klon,klev)
321	C Tendances Ondes de G non oro (runs strato).
322	REAL d_u_hin(klon,klev), d_v_hin(klon,klev)
323	REAL d_t_hin(klon,klev)
324
325	c Tendencies due to radiative scheme [K/s]
326	c d_t_rad,dtsw,dtlw,d_t_nirco2,d_t_nlte,d_t_euv
327	c are not computed at each physical timestep
328	c therefore, they are defined and saved in phys_state_var_mod
329
330	c Tendencies due to molecular viscosity and conduction
331	real d_t_conduc(klon,klev) ! [K/s]
332	real d_u_molvis(klon,klev) ! (m/s) /s
333	real d_v_molvis(klon,klev) ! (m/s) /s
334
335	c Tendencies due to molecular diffusion
336	real d_q_moldif(klon,klev,nqmax)
337
338	c
339	c Variables liees a l'ecriture de la bande histoire physique
340	c
341	INTEGER ecrit_mth
342	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
343	c
344	INTEGER ecrit_day
345	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
346	c
347	INTEGER ecrit_ins
348	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
349	c
350	integer itau_w ! pas de temps ecriture = itap + itau_phy
351
352	c Variables locales pour effectuer les appels en serie
353	c
354	REAL t_seri(klon,klev)
355	REAL u_seri(klon,klev), v_seri(klon,klev)
356	c
357	REAL :: tr_seri(klon,klev,nqmax)
358	REAL :: d_tr(klon,klev,nqmax)
359
360	c pour ioipsl
361	INTEGER nid_day, nid_mth, nid_ins
362	SAVE nid_day, nid_mth, nid_ins
363	INTEGER nhori, nvert, idayref
364	REAL zsto, zout, zsto1, zsto2, zero
365	parameter (zero=0.0e0)
366	real zjulian
367	save zjulian
368
369	CHARACTER*2 str2
370	character*20 modname
371	character*80 abort_message
372	logical ok_sync
373
374	character*30 nom_fichier
375	character*10 varname
376	character*40 vartitle
377	character*20 varunits
378	C Variables liees au bilan d'energie et d'enthalpi
379	REAL ztsol(klon)
380	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
381	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
382	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
383	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
384	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
385	REAL d_h_vcol_phy
386	REAL fs_bound, fq_bound
387	SAVE d_h_vcol_phy
388	REAL zero_v(klon),zero_v2(klon,klev)
389	CHARACTER*15 ztit
390	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
391	SAVE ip_ebil
392	DATA ip_ebil/2/
393	INTEGER if_ebil ! level for energy conserv. dignostics
394	SAVE if_ebil
395	c+jld ec_conser
396	REAL d_t_ec(klon,klev) ! tendance du a la conversion Ec -> E thermique
397	c-jld ec_conser
398
399	c TEST VENUS...
400	REAL mang(klon,klev) ! moment cinetique
401	REAL mangtot ! moment cinetique total
402
403	c cell_area for outputs in hist*
404	REAL cell_area_out(klon)
405	#ifdef MESOSCALE
406	REAL :: dt_dyn(klev)
407	#endif
408	c Declaration des constantes et des fonctions thermodynamiques
409	c
410	#include "YOMCST.h"
411
412	c======================================================================
413	c INITIALISATIONS
414	c================
415
416	modname = 'physiq'
417	ok_sync=.TRUE.
418
419	bilansmc = 0
420	ballons = 0
421	! NE FONCTIONNENT PAS ENCORE EN PARALLELE !!!
422	#ifndef MESOSCALE
423	if (is_parallel) then
424	bilansmc = 0
425	ballons = 0
426	endif
427	#endif
428	IF (if_ebil.ge.1) THEN
429	DO i=1,klon
430	zero_v(i)=0.
431	END DO
432	DO i=1,klon
433	DO j=1,klev
434	zero_v2(i,j)=0.
435	END DO
436	END DO
437	END IF
438
439	c PREMIER APPEL SEULEMENT
440	c========================
441	IF (debut) THEN
442	allocate(source(klon,nqmax))
443
444	#ifdef CPP_XIOS
445	! Initialize XIOS context
446	write(,) "physiq: call wxios_context_init"
447	CALL wxios_context_init
448	#endif
449
450	! The call to suphec is now done in iniphysiq_mod (interface)
451	! CALL suphec ! initialiser constantes et parametres phys.
452
453	IF (if_ebil.ge.1) d_h_vcol_phy=0.
454	c
455	c appel a la lecture du physiq.def
456	c
457	call conf_phys(ok_journe, ok_mensuel,
458	. ok_instan,
459	. if_ebil)
460
461	call phys_state_var_init(nqmax)
462	c
463	c Initialising Hedin model for upper atm
464	c (to be revised when coupled to chemistry) :
465	call conc_init
466
467	! initialise physics counter
468
469	itap = 0
470
471	#ifdef MESOSCALE
472	print*,'check pdtphys',pdtphys
473	PRINT*,'check phisfi ',pphis(1),pphis(klon)
474	PRINT*,'check geop',pphi(1,1),pphi(klon,klev)
475	PRINT*,'check radsol',radsol(1),radsol(klon)
476	print*,'check ppk',ppk(1,1),ppk(klon,klev)
477	print*,'check ftsoil',ftsoil(1,1),ftsoil(klon,nsoilmx)
478	print*,'check ftsol',ftsol(1),ftsol(klon)
479	print*, "check temp", t(1,1),t(klon,klev)
480	print*, "check pres",paprs(1,1),paprs(klon,klev),pplay(1,1),
481	. pplay(klon,klev)
482	print*, "check u", u(1,1),u(klon,klev)
483	print*, "check v", v(1,1),v(klon,klev)
484	print*,'check falbe',falbe(1),falbe(klon)
485	!nqtot=nqmax
486	!ALLOCATE(tname(nqtot))
487	!tname=noms
488	zmea=0.
489	zstd=0.
490	zsig=0.
491	zgam=0.
492	zthe=0.
493	dtime=pdtphys
494	#else
495	c
496	c Lecture startphy.nc :
497	c
498	CALL phyetat0 ("startphy.nc")
499	IF (.not.startphy_file) THEN
500	! Additionnal academic initializations
501	ftsol(:)=t(:,1) ! surface temperature as in first atm. layer
502	DO isoil=1, nsoilmx
503	! subsurface temperatures equal to surface temperature
504	ftsoil(:,isoil)=ftsol(:)
505	ENDDO
506	ENDIF
507	#endif
508
509	c dtime est defini dans tabcontrol.h et lu dans startphy
510	c pdtphys est calcule a partir des nouvelles conditions:
511	c Reinitialisation du pas de temps physique quand changement
512	IF (ABS(dtime-pdtphys).GT.0.001) THEN
513	WRITE(lunout,*) 'Pas physique a change',dtime,
514	. pdtphys
515	c abort_message='Pas physique n est pas correct '
516	c call abort_physic(modname,abort_message,1)
517	c----------------
518	c pour initialiser convenablement le time_counter, il faut tenir compte
519	c du changement de dtime en changeant itau_phy (point de depart)
520	itau_phy = NINT(itau_phy*dtime/pdtphys)
521	c----------------
522	dtime=pdtphys
523	ENDIF
524
525	radpas = NINT(RDAY/pdtphys/nbapp_rad)
526
527	CALL printflag( ok_journe,ok_instan )
528
529	#ifdef CPP_XIOS
530	write(,) "physiq: call initialize_xios_output"
531	call initialize_xios_output(rjourvrai,gmtime,pdtphys,RDAY,
532	& presnivs,pseudoalt)
533	#endif
534
535	c
536	c---------
537	c FLOTT
538	IF (ok_orodr) THEN
539	DO i=1,klon
540	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
541	ENDDO
542	CALL SUGWD(klon,klev,paprs,pplay)
543	DO i=1,klon
544	zuthe(i)=0.
545	zvthe(i)=0.
546	if(zstd(i).gt.10.)then
547	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
548	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
549	endif
550	ENDDO
551	ENDIF
552
553	if (bilansmc.eq.1) then
554	C OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
555	C DU BILAN DE MOMENT ANGULAIRE.
556	open(27,file='aaam_bud.out',form='formatted')
557	open(28,file='fields_2d.out',form='formatted')
558	write(,)'Ouverture de aaam_bud.out (FL Vous parle)'
559	write(,)'Ouverture de fields_2d.out (FL Vous parle)'
560	endif !bilansmc
561
562	c--------------SLEBONNOIS
563	C OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
564	C DES BALLONS
565	if (ballons.eq.1) then
566	open(30,file='ballons-lat.out',form='formatted')
567	open(31,file='ballons-lon.out',form='formatted')
568	open(32,file='ballons-u.out',form='formatted')
569	open(33,file='ballons-v.out',form='formatted')
570	open(34,file='ballons-alt.out',form='formatted')
571	write(,)'Ouverture des ballons*.out'
572	endif !ballons
573	c-------------
574
575	c---------
576	C TRACEURS
577	C source dans couche limite
578	source(:,:) = 0.0 ! pas de source, pour l'instant
579	c---------
580
581	c---------
582	c INITIALIZE THERMOSPHERIC PARAMETERS
583	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
584
585	if (callthermos) then
586	if(solvarmod.eq.0) call param_read
587	if(solvarmod.eq.1) call param_read_e107
588	endif
589
590	c Initialisation (recomputed in concentration2)
591	do ig=1,klon
592	do j=1,klev
593	rnew(ig,j)=R
594	cpnew(ig,j)=cpdet(t(ig,j))
595	mmean(ig,j)=RMD
596	akknew(ig,j)=1.e-4
597	rho(ig,j)=pplay(ig,j)mmean(ig,j)1e-3/(rnew(ig,j)*t(ig,j))
598	enddo
599	c stop
600
601	enddo
602
603	IF(callthermos.or.callnlte.or.callnirco2) THEN
604	call compo_hedin83_init2
605	ENDIF
606	if (callnlte.and.nltemodel.eq.2) call nlte_setup
607	if (callnirco2.and.nircorr.eq.1) call nir_leedat
608	c---------
609
610	c
611	c Verifications:
612	c
613	IF (nlon .NE. klon) THEN
614	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
615	. klon
616	abort_message='nlon et klon ne sont pas coherents'
617	call abort_physic(modname,abort_message,1)
618	ENDIF
619	IF (nlev .NE. klev) THEN
620	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
621	. klev
622	abort_message='nlev et klev ne sont pas coherents'
623	call abort_physic(modname,abort_message,1)
624	ENDIF
625	c
626	IF (dtimeREAL(radpas).GT.(RDAY0.25).AND.cycle_diurne)
627	$ THEN
628	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
629	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
630	abort_message='Nbre d appels au rayonnement insuffisant'
631	call abort_physic(modname,abort_message,1)
632	ENDIF
633	c
634	WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
635	. iflag_ajs
636	c
637	ecrit_mth = NINT(RDAY/dtime*ecriphy) ! tous les ecritphy jours
638	IF (ok_mensuel) THEN
639	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
640	. ecrit_mth
641	ENDIF
642
643	ecrit_day = NINT(RDAY/dtime *1.0) ! tous les jours
644	IF (ok_journe) THEN
645	WRITE(lunout,*)'La frequence de sortie journaliere est de ',
646	. ecrit_day
647	ENDIF
648
649	ecrit_ins = NINT(RDAY/dtime*ecriphy) ! Fraction de jour reglable
650	IF (ok_instan) THEN
651	WRITE(lunout,*)'La frequence de sortie instant. est de ',
652	. ecrit_ins
653	ENDIF
654
655	c Initialisation des sorties
656	c===========================
657
658	#ifdef CPP_IOIPSL
659
660	#ifdef histhf
661	#include "ini_histhf.h"
662	#endif
663
664	#ifdef histday
665	#include "ini_histday.h"
666	#endif
667
668	#ifdef histmth
669	#include "ini_histmth.h"
670	#endif
671
672	#ifdef histins
673	#include "ini_histins.h"
674	#endif
675
676	#endif
677
678	c
679	c Initialiser les valeurs de u pour calculs tendances
680	c (pour T, c'est fait dans phyetat0)
681	c
682	DO k = 1, klev
683	DO i = 1, klon
684	u_ancien(i,k) = u(i,k)
685	ENDDO
686	ENDDO
687
688	c---------
689	c Ecriture fichier initialisation
690	c PRINT*,'Ecriture Initial_State.csv'
691	c OPEN(88,file='Trac_Point.csv',
692	c & form='formatted')
693	c---------
694
695	c---------
696	c Initialisation des parametres des nuages
697	c===============================================
698
699	c MICROPHY SANS CHIMIE: seulement si full microphy (cl_scheme=2)
700
701	if (ok_chem.and..not.ok_cloud) then
702	print*,"LA CHIMIE A BESOIN DE LA MICROPHYSIQUE"
703	print*,"ok_cloud doit etre = a ok_chem"
704	stop
705	endif
706	if (.not.ok_chem.and.ok_cloud.and.(cl_scheme.eq.1)) then
707	print*,"cl_scheme=1 doesnot work without chemistry"
708	stop
709	endif
710	if (.not.ok_chem.and.ok_cloud.and.(cl_scheme.eq.2)) then
711	print*,"Full microphysics without chemistry"
712	c indexation of microphys tracers
713	CALL chemparam_ini()
714	endif
715
716	! number of microphysical tracers
717
718	nmicro = 0
719	if (ok_cloud .and. (cl_scheme == 1)) nmicro = 2
720	if (ok_cloud .and. (cl_scheme == 2)) nmicro = 12
721
722	c CAS 1D POUR MICROPHYS Aurelien
723	if ((nlon .EQ. 1) .AND. ok_cloud .AND. (cl_scheme.eq.1)) then
724	PRINT*,'Open profile_cloud_parameters.csv'
725	OPEN(66,file='profile_cloud_parameters.csv',
726	& form='formatted')
727	endif
728
729	if ((nlon .EQ. 1) .AND. ok_sedim .AND. (cl_scheme.eq.1)) then
730	PRINT*,'Open profile_cloud_sedim.csv'
731	OPEN(77,file='profile_cloud_sedim.csv',
732	& form='formatted')
733	endif
734
735	c INIT PHOTOCHEMISTRY ! includes the indexation of microphys tracers
736	c if ((nlon .GT. 1) .AND. ok_chem) then
737	c !!! DONC 3D !!! POURQUOI ???
738	if (ok_chem) then
739	CALL chemparam_ini()
740	endif
741
742	c INIT MICROPHYS SCHEME 1 (AURELIEN)
743	if ((nlon .GT. 1) .AND. ok_cloud .AND. (cl_scheme.eq.1)) then
744	c !!! DONC 3D !!!
745	CALL cloud_ini(nlon,nlev)
746	endif
747
748	ENDIF ! debut
749
750	! ------------------------------------------------------
751	! Initializations done at every physical timestep:
752	! ------------------------------------------------------
753
754	c Mettre a zero des variables de sortie (pour securite)
755	c
756	DO i = 1, klon
757	d_ps(i) = 0.0
758	ENDDO
759	DO k = 1, klev
760	DO i = 1, klon
761	d_t(i,k) = 0.0
762	d_u(i,k) = 0.0
763	d_v(i,k) = 0.0
764	ENDDO
765	ENDDO
766	DO iq = 1, nqmax
767	DO k = 1, klev
768	DO i = 1, klon
769	d_qx(i,k,iq) = 0.0
770	ENDDO
771	ENDDO
772	ENDDO
773	c
774	c Ne pas affecter les valeurs entrees de u, v, h, et q
775	c
776	DO k = 1, klev
777	DO i = 1, klon
778	t_seri(i,k) = t(i,k)
779	u_seri(i,k) = u(i,k)
780	v_seri(i,k) = v(i,k)
781	ENDDO
782	ENDDO
783	DO iq = 1, nqmax
784	DO k = 1, klev
785	DO i = 1, klon
786	tr_seri(i,k,iq) = qx(i,k,iq)
787	ENDDO
788	ENDDO
789	ENDDO
790	C
791	DO i = 1, klon
792	ztsol(i) = ftsol(i)
793	ENDDO
794	C
795	IF (if_ebil.ge.1) THEN
796	ztit='after dynamic'
797	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
798	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
799	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
800	C Comme les tendances de la physique sont ajoute dans la dynamique,
801	C on devrait avoir que la variation d'entalpie par la dynamique
802	C est egale a la variation de la physique au pas de temps precedent.
803	C Donc la somme de ces 2 variations devrait etre nulle.
804	call diagphy(cell_area,ztit,ip_ebil
805	e , zero_v, zero_v, zero_v, zero_v, zero_v
806	e , zero_v, zero_v, zero_v, ztsol
807	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
808	s , fs_bound, fq_bound )
809	END IF
810
811	c====================================================================
812	c XIOS outputs
813
814	#ifdef CPP_XIOS
815	! update XIOS time/calendar
816	call update_xios_timestep
817	#endif
818
819	c====================================================================
820	c Diagnostiquer la tendance dynamique
821	c
822	IF (ancien_ok) THEN
823	DO k = 1, klev
824	DO i = 1, klon
825	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
826	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
827	ENDDO
828	ENDDO
829
830	! ADAPTATION GCM POUR CP(T)
831	do i=1,klon
832	flux_dyn(i,1) = 0.0
833	do j=2,klev
834	flux_dyn(i,j) = flux_dyn(i,j-1)
835	. +cpdet(t_seri(i,j-1))/RGd_t_dyn(i,j-1)(paprs(i,j-1)-paprs(i,j))
836	enddo
837	enddo
838
839	ELSE
840	DO k = 1, klev
841	DO i = 1, klon
842	d_u_dyn(i,k) = 0.0
843	d_t_dyn(i,k) = 0.0
844	ENDDO
845	ENDDO
846	ancien_ok = .TRUE.
847	ENDIF
848	c====================================================================
849
850	c Calcule de vitesse verticale a partir de flux de masse verticale
851	DO k = 1, klev
852	DO i = 1, klon
853	omega(i,k) = RG*flxmw(i,k) / cell_area(i)
854	END DO
855	END DO
856
857	c======
858	c GEOP CORRECTION
859	c
860	c Ajouter le geopotentiel du sol:
861	c
862	DO k = 1, klev
863	DO i = 1, klon
864	zphi(i,k) = pphi(i,k) + pphis(i)
865	ENDDO
866	ENDDO
867
868	c............................
869	c CETTE CORRECTION VA DE PAIR AVEC DES MODIFS DE LEAPFROG(_p)
870	c ELLE MARCHE A 50 NIVEAUX (si mmean constante...)
871	c MAIS PAS A 78 NIVEAUX (quand mmean varie...)
872	c A ANALYSER PLUS EN DETAIL AVANT D'UTILISER
873	c............................
874	c zphi is recomputed (pphi is not ok if mean molecular mass varies)
875	c with dphi = RT/mmean d(ln p) [evaluated at interface]
876
877	c DO i = 1, klon
878	c zphi(i,1) = pphis(i) + Rt_seri(i,1)/mmean(i,1)1000.
879	c * *( log(paprs(i,1)) - log(pplay(i,1)) )
880	c DO k = 2, klev
881	c zphi(i,k) = zphi(i,k-1)
882	c * + R500.(t_seri(i,k)/mmean(i,k)+t_seri(i,k-1)/mmean(i,k-1))
883	c * * (log(pplay(i,k-1)) - log(pplay(i,k)))
884	c ENDDO
885	c ENDDO
886	c............................
887	c=====
888
889	c calcul du geopotentiel aux niveaux intercouches
890	c ponderation des altitudes au niveau des couches en dp/p
891
892	DO k=1,klev
893	DO i=1,klon
894	zzlay(i,k)=zphi(i,k)/RG ! [m]
895	ENDDO
896	ENDDO
897	DO i=1,klon
898	zzlev(i,1)=pphis(i)/RG ! [m]
899	ENDDO
900	DO k=2,klev
901	DO i=1,klon
902	z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
903	z2=(paprs(i,k) +pplay(i,k))/(paprs(i,k) -pplay(i,k))
904	zzlev(i,k)=(z1zzlay(i,k-1)+z2zzlay(i,k))/(z1+z2)
905	ENDDO
906	ENDDO
907	DO i=1,klon
908	zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
909	ENDDO
910
911	c====================================================================
912	c
913	c Verifier les temperatures
914	c
915	CALL hgardfou(t_seri,ftsol,'debutphy')
916
917	c====================================================================
918	c Orbite et eclairement
919	c====================================================================
920
921	c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0.
922	c donc pas de variations de Ls, ni de dist.
923	c La seule chose qui compte, c'est la rotation de la planete devant
924	c le Soleil...
925
926	zlongi = 0.0
927	dist = 0.72 ! en UA
928
929	c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite
930	c a sa valeur, et prendre en compte leur evolution,
931	c il faudra refaire un orbite.F...
932	c CALL orbite(zlongi,dist)
933
934	IF (cycle_diurne) THEN
935	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
936	CALL zenang(zlongi,gmtime,zdtime,latitude_deg,longitude_deg,
937	& rmu0,fract)
938	ELSE
939	call mucorr(klon,zlongi,latitude_deg,rmu0,fract)
940	ENDIF
941
942	! print fraction of venus day
943
944	if (is_master) then
945	print*, 'gmtime = ', gmtime
946	end if
947
948	if (iflag_trac.eq.1) then
949	!====================================================================
950	! Case 1: pseudo-chemistry with relaxation toward fixed profile
951	!====================================================================
952	if (tr_scheme.eq.1) then
953
954	call phytrac_relax (debut,lafin,nqmax,
955	I nlon,nlev,dtime,pplay,
956	O tr_seri)
957
958	elseif (tr_scheme.eq.2) then
959	!====================================================================
960	! Case 2: surface emission
961	! For the moment, inspired from Mars version
962	! However, the variable 'source' could be used in physiq
963	! so the call to phytrac_emiss could be to initialise it.
964	!====================================================================
965
966	call phytrac_emiss ( (rjourvrai+gmtime)*RDAY,
967	I debut,lafin,nqmax,
968	I nlon,nlev,dtime,paprs,
969	I latitude_deg,longitude_deg,
970	O tr_seri)
971
972	else if (tr_scheme.eq.3) then
973	!====================================================================
974	! Case 3: Full chemistry and/or clouds.
975	! routines are called every "chempas" physical timestep.
976	! if the physics is called 96000 times per venus day,
977	! then chempas = 4
978	!====================================================================
979
980	nbapp_chem = 24000 ! corresponds to 420 seconds
981	chempas = nint(rday/pdtphys/nbapp_chem)
982	zctime = dtime*real(chempas) ! chemical timestep (420s)
983
984	if (mod(itap,chempas) == 0) then ! <------- start of chemistry supercycling
985
986	! photochemistry and microphysics
987
988	call phytrac_chimie(debut,
989	$ gmtime,
990	$ nqmax,
991	$ klon,
992	$ latitude_deg,
993	$ longitude_deg,
994	$ nlev,
995	$ zctime,
996	$ t_seri,
997	$ pplay,
998	$ tr_seri,
999	$ d_tr_chem,
1000	$ iter)
1001
1002	if (ok_sedim) then
1003	if (cl_scheme == 1) then
1004
1005	! sedimentation for simplified microphysics
1006
1007	#ifndef MESOSCALE
1008	call new_cloud_sedim(klon,
1009	$ nlev,
1010	$ zctime,
1011	$ pplay,
1012	$ paprs,
1013	$ t_seri,
1014	$ tr_seri,
1015	$ d_tr_chem,
1016	$ d_tr_sed(:,:,1:2),
1017	$ d_tr_ssed,
1018	$ nqmax,
1019	$ Fsedim)
1020
1021	! test to avoid nans
1022
1023	do k = 1, klev
1024	do i = 1, klon
1025	if ((d_tr_sed(i,k,1) /= d_tr_sed(i,k,1)) .or.
1026	$ (d_tr_sed(i,k,2) /= d_tr_sed(i,k,2))) then
1027	print*,'sedim NaN PROBLEM'
1028	print*,'d_tr_sed Nan?',d_tr_sed(i,k,:)
1029	print*,'Temp',t_seri(i,k)
1030	print*,'lat-lon',i,'level',k,'zctime',zctime
1031	print*,'F_sed',Fsedim(i,k)
1032	d_tr_sed(i,k,:) = 0.
1033	end if
1034	end do
1035	end do
1036
1037	! tendency due to sedimentation
1038
1039	d_tr_sed(:,:,:) = d_tr_sed(:,:,:)/zctime
1040	Fsedim(:,1:klev) = Fsedim(:,1:klev)/zctime
1041	Fsedim(:,klev+1) = 0.
1042
1043	else if (cl_scheme == 2) then
1044
1045	! sedimentation for detailed microphysics
1046
1047	d_tr_sed(:,:,:) = 0.
1048
1049	do i = 1, klon
1050
1051	! mode 1
1052
1053	m0_mode1(:,1) = tr_seri(i,:,i_m0_mode1drop)
1054	m0_mode1(:,2) = tr_seri(i,:,i_m0_mode1ccn)
1055	m3_mode1(:,1) = tr_seri(i,:,i_m3_mode1sa)
1056	m3_mode1(:,2) = tr_seri(i,:,i_m3_mode1w)
1057	m3_mode1(:,3) = tr_seri(i,:,i_m3_mode1ccn)
1058
1059	call drop_sedimentation(zctime,klev,m0_mode1,m3_mode1,
1060	$ paprs(i,:),zzlev(i,:),
1061	$ zzlay(i,:),t_seri(i,:),1,
1062	$ d_ccn_sed(:,1),d_drop_sed,
1063	$ d_ccn_sed(:,2),d_liq_sed)
1064
1065	d_tr_sed(i,:,i_m0_mode1drop)= d_tr_sed(i,:,i_m0_mode1drop)
1066	$ + d_drop_sed
1067	d_tr_sed(i,:,i_m0_mode1ccn) = d_tr_sed(i,:,i_m0_mode1ccn)
1068	$ + d_ccn_sed(:,1)
1069	d_tr_sed(i,:,i_m3_mode1ccn) = d_tr_sed(i,:,i_m3_mode1ccn)
1070	$ + d_ccn_sed(:,2)
1071	d_tr_sed(i,:,i_m3_mode1sa) = d_tr_sed(i,:,i_m3_mode1sa)
1072	$ + d_liq_sed(:,1)
1073	d_tr_sed(i,:,i_m3_mode1w) = d_tr_sed(i,:,i_m3_mode1w)
1074	$ + d_liq_sed(:,2)
1075
1076	! mode 2
1077
1078	m0_mode2(:,1) = tr_seri(i,:,i_m0_mode2drop)
1079	m0_mode2(:,2) = tr_seri(i,:,i_m0_mode2ccn)
1080	m3_mode2(:,1) = tr_seri(i,:,i_m3_mode2sa)
1081	m3_mode2(:,2) = tr_seri(i,:,i_m3_mode2w)
1082	m3_mode2(:,3) = tr_seri(i,:,i_m3_mode2ccn)
1083
1084	call drop_sedimentation(zctime,klev,m0_mode2,m3_mode2,
1085	$ paprs(i,:),zzlev(i,:),
1086	$ zzlay(i,:),t_seri(i,:),2,
1087	$ d_ccn_sed(:,1),d_drop_sed,
1088	$ d_ccn_sed(:,2),d_liq_sed)
1089
1090	d_tr_sed(i,:,i_m0_mode2drop)= d_tr_sed(i,:,i_m0_mode2drop)
1091	$ + d_drop_sed
1092	d_tr_sed(i,:,i_m0_mode2ccn) = d_tr_sed(i,:,i_m0_mode2ccn)
1093	$ + d_ccn_sed(:,1)
1094	d_tr_sed(i,:,i_m3_mode2ccn) = d_tr_sed(i,:,i_m3_mode2ccn)
1095	$ + d_ccn_sed(:,2)
1096	d_tr_sed(i,:,i_m3_mode2sa) = d_tr_sed(i,:,i_m3_mode2sa)
1097	$ + d_liq_sed(:,1)
1098	d_tr_sed(i,:,i_m3_mode2w) = d_tr_sed(i,:,i_m3_mode2w)
1099	$ + d_liq_sed(:,2)
1100
1101	! aer
1102
1103	call aer_sedimentation(zctime,klev,
1104	$ tr_seri(i,:,i_m0_aer),
1105	$ tr_seri(i,:,i_m3_aer),
1106	$ paprs(i,:),zzlev(i,:),
1107	$ zzlay(i,:),t_seri(i,:),
1108	$ d_tr_sed(i,:,i_m0_aer),
1109	$ d_tr_sed(i,:,i_m3_aer),
1110	$ aer_flux)
1111
1112	end do
1113
1114	! tendency due to sedimentation
1115
1116	do iq = nqmax-nmicro+1,nqmax
1117	d_tr_sed(:,:,iq) = d_tr_sed(:,:,iq)/zctime
1118	end do
1119	#endif
1120	end if ! cl_scheme
1121	end if ! ok_sedim
1122	end if ! mod(itap,chempas) <------- end of chemistry supercycling
1123
1124	! update tracers (chemistry)
1125
1126	do iq = 1, nqmax - nmicro
1127	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_chem(:,:,iq)*dtime
1128	end do
1129
1130	! update tracers (sedimentation)
1131
1132	if (ok_sedim) then
1133	if (cl_scheme == 1) then
1134	tr_seri(:,:,i_h2so4liq) = tr_seri(:,:,i_h2so4liq)
1135	$ + d_tr_sed(:,:,1)*dtime
1136	tr_seri(:,:,i_h2oliq) = tr_seri(:,:,i_h2oliq)
1137	$ + d_tr_sed(:,:,2)*dtime
1138	else if (cl_scheme == 2) then
1139	do iq = nqmax-nmicro+1,nqmax
1140	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_sed(:,:,iq)*dtime
1141	end do
1142	end if ! cl_scheme
1143	end if ! ok_sedim
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