Context Navigation

physiq_mod.F @ 2135

Last change on this file since 2135 was 2135, checked in by slebonnois, 6 years ago
SL, Venus: new keys for flexibility cp0/cp(T) and Held-Suarez type physics
File size: 63.7 KB

Line
1	!
2	! $Id: $
3	!
4	MODULE physiq_mod
5
6	IMPLICIT NONE
7
8	CONTAINS
9
10	SUBROUTINE physiq (nlon,nlev,nqmax,
11	. debut,lafin,rjourvrai,gmtime,pdtphys,
12	. paprs,pplay,ppk,pphi,pphis,presnivs,
13	. u,v,t,qx,
14	. flxmw,
15	. d_u, d_v, d_t, d_qx, d_ps)
16
17	c======================================================================
18	c
19	c Modifications pour la physique de Venus
20	c S. Lebonnois (LMD/CNRS) Septembre 2005
21	c
22	c ---------------------------------------------------------------------
23	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
24	c
25	c Objet: Moniteur general de la physique du modele
26	cAA Modifications quant aux traceurs :
27	cAA - uniformisation des parametrisations ds phytrac
28	cAA - stockage des moyennes des champs necessaires
29	cAA en mode traceur off-line
30	c modif ( P. Le Van , 12/10/98 )
31	c
32	c Arguments:
33	c
34	c nlon----input-I-nombre de points horizontaux
35	c nlev----input-I-nombre de couches verticales
36	c nqmax---input-I-nombre de traceurs
37	c debut---input-L-variable logique indiquant le premier passage
38	c lafin---input-L-variable logique indiquant le dernier passage
39	c rjour---input-R-numero du jour de l'experience
40	c gmtime--input-R-fraction de la journee (0 a 1)
41	c pdtphys-input-R-pas d'integration pour la physique (seconde)
42	c paprs---input-R-pression pour chaque inter-couche (en Pa)
43	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
44	c ppk ---input-R-fonction d'Exner au milieu de couche
45	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
46	c pphis---input-R-geopotentiel du sol
47	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
48	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
49	c v-------input-R-vitesse Y (de S a N) en m/s
50	c t-------input-R-temperature (K)
51	c qx------input-R-mass mixing ratio traceurs (kg/kg)
52	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
53	c flxmw---input-R-flux de masse vertical en kg/s
54	c
55	c d_u-----output-R-tendance physique de "u" (m/s/s)
56	c d_v-----output-R-tendance physique de "v" (m/s/s)
57	c d_t-----output-R-tendance physique de "t" (K/s)
58	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
59	c d_ps----output-R-tendance physique de la pression au sol
60	c======================================================================
61	USE ioipsl
62	! USE histcom ! not needed; histcom is included in ioipsl
63	use dimphy
64	USE geometry_mod,only: longitude, latitude, ! in radians
65	& longitude_deg,latitude_deg, ! in degrees
66	& cell_area,dx,dy
67	USE phys_state_var_mod ! Variables sauvegardees de la physique
68	USE cpdet_phy_mod, only: cpdet, t2tpot
69	USE chemparam_mod
70	USE conc
71	USE 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	#endif
94	IMPLICIT none
95	c======================================================================
96	c CLEFS CPP POUR LES IO
97	c =====================
98	#ifndef MESOSCALE
99	c#define histhf
100	#define histday
101	#define histmth
102	#define histins
103	#endif
104	c======================================================================
105	#include "dimsoil.h"
106	#include "clesphys.h"
107	#include "iniprint.h"
108	#include "timerad.h"
109	#include "tabcontrol.h"
110	#include "nirdata.h"
111	#include "hedin.h"
112	c======================================================================
113	LOGICAL ok_journe ! sortir le fichier journalier
114	save ok_journe
115	c PARAMETER (ok_journe=.true.)
116	c
117	LOGICAL ok_mensuel ! sortir le fichier mensuel
118	save ok_mensuel
119	c PARAMETER (ok_mensuel=.true.)
120	c
121	LOGICAL ok_instan ! sortir le fichier instantane
122	save ok_instan
123	c PARAMETER (ok_instan=.true.)
124	c
125	c======================================================================
126	c
127	c Variables argument:
128	c
129	INTEGER nlon
130	INTEGER nlev
131	INTEGER nqmax
132	REAL rjourvrai
133	REAL gmtime
134	REAL pdtphys
135	LOGICAL debut, lafin
136	REAL paprs(klon,klev+1)
137	REAL pplay(klon,klev)
138	REAL pphi(klon,klev)
139	REAL pphis(klon)
140	REAL presnivs(klev)
141
142	! ADAPTATION GCM POUR CP(T)
143	REAL ppk(klon,klev)
144
145	REAL u(klon,klev)
146	REAL v(klon,klev)
147	REAL t(klon,klev)
148	REAL qx(klon,klev,nqmax)
149
150	REAL d_u_dyn(klon,klev)
151	REAL d_t_dyn(klon,klev)
152
153	REAL flxmw(klon,klev)
154
155	REAL d_u(klon,klev)
156	REAL d_v(klon,klev)
157	REAL d_t(klon,klev)
158	REAL d_qx(klon,klev,nqmax)
159	REAL d_ps(klon)
160
161	logical ok_hf
162	real ecrit_hf
163	integer nid_hf
164	save ok_hf, ecrit_hf, nid_hf
165
166	#ifdef histhf
167	data ok_hf,ecrit_hf/.true.,0.25/
168	#else
169	data ok_hf/.false./
170	#endif
171
172	c Variables propres a la physique
173	c
174	INTEGER,save :: itap ! compteur pour la physique
175	REAL delp(klon,klev) ! epaisseur d'une couche
176	REAL omega(klon,klev) ! vitesse verticale en Pa/s
177
178
179	INTEGER igwd,idx(klon),itest(klon)
180	c
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
278	INTEGER itaprad
279	SAVE itaprad
280	c
281	REAL dist, rmu0(klon), fract(klon)
282	REAL zdtime, zlongi
283	c
284	INTEGER i, k, iq, ig, j, ll, ilon, ilat, ilev, isoil
285	c
286	REAL zphi(klon,klev)
287	REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
288	real tsurf(klon)
289
290	c va avec nlte_tcool
291	INTEGER ierr_nlte
292	REAL varerr
293
294	c Variables du changement
295	c
296	c ajs: ajustement sec
297	c vdf: couche limite (Vertical DiFfusion)
298	REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax)
299	REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
300	c
301	REAL d_ts(klon)
302	c
303	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
304	REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax)
305	c
306	CMOD LOTT: Tendances Orography Sous-maille
307	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
308	REAL d_t_oro(klon,klev)
309	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
310	REAL d_t_lif(klon,klev)
311	C Tendances Ondes de G non oro (runs strato).
312	REAL d_u_hin(klon,klev), d_v_hin(klon,klev)
313	REAL d_t_hin(klon,klev)
314
315	c Tendencies due to radiative scheme [K/s]
316	c d_t_rad,dtsw,dtlw,d_t_nirco2,d_t_nlte,d_t_euv
317	c are not computed at each physical timestep
318	c therefore, they are defined and saved in phys_state_var_mod
319
320	c Tendencies due to molecular viscosity and conduction
321	real d_t_conduc(klon,klev) ! [K/s]
322	real d_u_molvis(klon,klev) ! (m/s) /s
323	real d_v_molvis(klon,klev) ! (m/s) /s
324
325	c Tendencies due to molecular diffusion
326	real d_q_moldif(klon,klev,nqmax)
327
328	c
329	c Variables liees a l'ecriture de la bande histoire physique
330	c
331	INTEGER ecrit_mth
332	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
333	c
334	INTEGER ecrit_day
335	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
336	c
337	INTEGER ecrit_ins
338	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
339	c
340	integer itau_w ! pas de temps ecriture = itap + itau_phy
341
342	c Variables locales pour effectuer les appels en serie
343	c
344	REAL t_seri(klon,klev)
345	REAL u_seri(klon,klev), v_seri(klon,klev)
346	c
347	REAL :: tr_seri(klon,klev,nqmax)
348	REAL :: d_tr(klon,klev,nqmax)
349
350	c Variables tendance sedimentation
351
352	REAL :: m0_mode1(klev,2),m0_mode2(klev,2)
353	REAL :: m3_mode1(klev,3),m3_mode2 (klev,3)
354	REAL :: d_drop_sed(klev),d_ccn_sed(klev,2),d_liq_sed(klev,2)
355	REAL :: aer_flux(klev)
356	REAL :: d_tr_sed(klon,klev,nqmax)
357	REAL :: d_tr_ssed(klon)
358	c
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
461	c
462	c Initialising Hedin model for upper atm
463	c (to be revised when coupled to chemistry) :
464	call conc_init
465	c
466	c Initialiser les compteurs:
467	c
468	itap = 0
469	itaprad = 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
531	write(,) "physiq: call initialize_xios_output"
532	call initialize_xios_output(rjourvrai,gmtime,pdtphys,RDAY,
533	& presnivs,pseudoalt)
534	#endif
535
536	c
537	c---------
538	c FLOTT
539	IF (ok_orodr) THEN
540	DO i=1,klon
541	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
542	ENDDO
543	CALL SUGWD(klon,klev,paprs,pplay)
544	DO i=1,klon
545	zuthe(i)=0.
546	zvthe(i)=0.
547	if(zstd(i).gt.10.)then
548	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
549	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
550	endif
551	ENDDO
552	ENDIF
553
554	if (bilansmc.eq.1) then
555	C OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
556	C DU BILAN DE MOMENT ANGULAIRE.
557	open(27,file='aaam_bud.out',form='formatted')
558	open(28,file='fields_2d.out',form='formatted')
559	write(,)'Ouverture de aaam_bud.out (FL Vous parle)'
560	write(,)'Ouverture de fields_2d.out (FL Vous parle)'
561	endif !bilansmc
562
563	c--------------SLEBONNOIS
564	C OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
565	C DES BALLONS
566	if (ballons.eq.1) then
567	open(30,file='ballons-lat.out',form='formatted')
568	open(31,file='ballons-lon.out',form='formatted')
569	open(32,file='ballons-u.out',form='formatted')
570	open(33,file='ballons-v.out',form='formatted')
571	open(34,file='ballons-alt.out',form='formatted')
572	write(,)'Ouverture des ballons*.out'
573	endif !ballons
574	c-------------
575
576	c---------
577	C TRACEURS
578	C source dans couche limite
579	source(:,:) = 0.0 ! pas de source, pour l'instant
580	c---------
581
582	c---------
583	c INITIALIZE THERMOSPHERIC PARAMETERS
584	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
585
586	if (callthermos) then
587	if(solvarmod.eq.0) call param_read
588	if(solvarmod.eq.1) call param_read_e107
589	endif
590
591	c Initialisation (recomputed in concentration2)
592	do ig=1,klon
593	do j=1,klev
594	rnew(ig,j)=R
595	cpnew(ig,j)=cpdet(t(ig,j))
596	mmean(ig,j)=RMD
597	akknew(ig,j)=1.e-4
598	rho(ig,j)=pplay(ig,j)mmean(ig,j)1e-3/(rnew(ig,j)*t(ig,j))
599	enddo
600	c stop
601
602	enddo
603
604	IF(callthermos.or.callnlte.or.callnirco2) THEN
605	call compo_hedin83_init2
606	ENDIF
607	if (callnlte.and.nltemodel.eq.2) call nlte_setup
608	if (callnirco2.and.nircorr.eq.1) call nir_leedat
609	c---------
610
611	c
612	c Verifications:
613	c
614	IF (nlon .NE. klon) THEN
615	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
616	. klon
617	abort_message='nlon et klon ne sont pas coherents'
618	call abort_physic(modname,abort_message,1)
619	ENDIF
620	IF (nlev .NE. klev) THEN
621	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
622	. klev
623	abort_message='nlev et klev ne sont pas coherents'
624	call abort_physic(modname,abort_message,1)
625	ENDIF
626	c
627	IF (dtimeREAL(radpas).GT.(RDAY0.25).AND.cycle_diurne)
628	$ THEN
629	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
630	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
631	abort_message='Nbre d appels au rayonnement insuffisant'
632	call abort_physic(modname,abort_message,1)
633	ENDIF
634	c
635	WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
636	. iflag_ajs
637	c
638	ecrit_mth = NINT(RDAY/dtime*ecriphy) ! tous les ecritphy jours
639	IF (ok_mensuel) THEN
640	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
641	. ecrit_mth
642	ENDIF
643
644	ecrit_day = NINT(RDAY/dtime *1.0) ! tous les jours
645	IF (ok_journe) THEN
646	WRITE(lunout,*)'La frequence de sortie journaliere est de ',
647	. ecrit_day
648	ENDIF
649
650	ecrit_ins = NINT(RDAY/dtime*ecriphy) ! Fraction de jour reglable
651	IF (ok_instan) THEN
652	WRITE(lunout,*)'La frequence de sortie instant. est de ',
653	. ecrit_ins
654	ENDIF
655
656	c Initialisation des sorties
657	c===========================
658
659	#ifdef CPP_IOIPSL
660
661	#ifdef histhf
662	#include "ini_histhf.h"
663	#endif
664
665	#ifdef histday
666	#include "ini_histday.h"
667	#endif
668
669	#ifdef histmth
670	#include "ini_histmth.h"
671	#endif
672
673	#ifdef histins
674	#include "ini_histins.h"
675	#endif
676
677	#endif
678
679	c
680	c Initialiser les valeurs de u pour calculs tendances
681	c (pour T, c'est fait dans phyetat0)
682	c
683	DO k = 1, klev
684	DO i = 1, klon
685	u_ancien(i,k) = u(i,k)
686	ENDDO
687	ENDDO
688
689	c---------
690	c Ecriture fichier initialisation
691	c PRINT*,'Ecriture Initial_State.csv'
692	c OPEN(88,file='Trac_Point.csv',
693	c & form='formatted')
694	c---------
695
696	c---------
697	c Initialisation des parametres des nuages
698	c===============================================
699
700	c MICROPHY SANS CHIMIE: seulement si full microphy (cl_scheme=2)
701
702	if (ok_chem.and..not.ok_cloud) then
703	print*,"LA CHIMIE A BESOIN DE LA MICROPHYSIQUE"
704	print*,"ok_cloud doit etre = a ok_chem"
705	stop
706	endif
707	if (.not.ok_chem.and.ok_cloud.and.(cl_scheme.eq.1)) then
708	print*,"cl_scheme=1 doesnot work without chemistry"
709	stop
710	endif
711	if (.not.ok_chem.and.ok_cloud.and.(cl_scheme.eq.2)) then
712	print*,"Full microphysics without chemistry"
713	c indexation of microphys tracers
714	CALL chemparam_ini()
715	endif
716
717	c number of microphysical tracers
718	nmicro=0
719	if (ok_cloud .AND. (cl_scheme.eq.1)) nmicro=2
720	if (ok_cloud .AND. (cl_scheme.eq.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	c======================================================================
750	c======================================================================
751	! ------------------------------------------------------
752	! Initializations done at every physical timestep:
753	! ------------------------------------------------------
754
755	c Mettre a zero des variables de sortie (pour securite)
756	c
757	DO i = 1, klon
758	d_ps(i) = 0.0
759	ENDDO
760	DO k = 1, klev
761	DO i = 1, klon
762	d_t(i,k) = 0.0
763	d_u(i,k) = 0.0
764	d_v(i,k) = 0.0
765	ENDDO
766	ENDDO
767	DO iq = 1, nqmax
768	DO k = 1, klev
769	DO i = 1, klon
770	d_qx(i,k,iq) = 0.0
771	ENDDO
772	ENDDO
773	ENDDO
774	c
775	c Ne pas affecter les valeurs entrees de u, v, h, et q
776	c
777	DO k = 1, klev
778	DO i = 1, klon
779	t_seri(i,k) = t(i,k)
780	u_seri(i,k) = u(i,k)
781	v_seri(i,k) = v(i,k)
782	ENDDO
783	ENDDO
784	DO iq = 1, nqmax
785	DO k = 1, klev
786	DO i = 1, klon
787	tr_seri(i,k,iq) = qx(i,k,iq)
788	ENDDO
789	ENDDO
790	ENDDO
791	C
792	DO i = 1, klon
793	ztsol(i) = ftsol(i)
794	ENDDO
795	C
796	IF (if_ebil.ge.1) THEN
797	ztit='after dynamic'
798	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
799	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
800	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
801	C Comme les tendances de la physique sont ajoute dans la dynamique,
802	C on devrait avoir que la variation d'entalpie par la dynamique
803	C est egale a la variation de la physique au pas de temps precedent.
804	C Donc la somme de ces 2 variations devrait etre nulle.
805	call diagphy(cell_area,ztit,ip_ebil
806	e , zero_v, zero_v, zero_v, zero_v, zero_v
807	e , zero_v, zero_v, zero_v, ztsol
808	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
809	s , fs_bound, fq_bound )
810	END IF
811
812	c====================================================================
813	c XIOS outputs
814
815	#ifdef CPP_XIOS
816	! update XIOS time/calendar
817	call update_xios_timestep
818	#endif
819
820	c====================================================================
821	c Diagnostiquer la tendance dynamique
822	c
823	IF (ancien_ok) THEN
824	DO k = 1, klev
825	DO i = 1, klon
826	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
827	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
828	ENDDO
829	ENDDO
830
831	! ADAPTATION GCM POUR CP(T)
832	do i=1,klon
833	flux_dyn(i,1) = 0.0
834	do j=2,klev
835	flux_dyn(i,j) = flux_dyn(i,j-1)
836	. +cpdet(t_seri(i,j-1))/RGd_t_dyn(i,j-1)(paprs(i,j-1)-paprs(i,j))
837	enddo
838	enddo
839
840	ELSE
841	DO k = 1, klev
842	DO i = 1, klon
843	d_u_dyn(i,k) = 0.0
844	d_t_dyn(i,k) = 0.0
845	ENDDO
846	ENDDO
847	ancien_ok = .TRUE.
848	ENDIF
849	c====================================================================
850
851	c Calcule de vitesse verticale a partir de flux de masse verticale
852	DO k = 1, klev
853	DO i = 1, klon
854	omega(i,k) = RG*flxmw(i,k) / cell_area(i)
855	END DO
856	END DO
857
858	c======
859	c GEOP CORRECTION
860	c
861	c Ajouter le geopotentiel du sol:
862	c
863	DO k = 1, klev
864	DO i = 1, klon
865	zphi(i,k) = pphi(i,k) + pphis(i)
866	ENDDO
867	ENDDO
868
869	c............................
870	c CETTE CORRECTION VA DE PAIR AVEC DES MODIFS DE LEAPFROG(_p)
871	c ELLE MARCHE A 50 NIVEAUX (si mmean constante...)
872	c MAIS PAS A 78 NIVEAUX (quand mmean varie...)
873	c A ANALYSER PLUS EN DETAIL AVANT D'UTILISER
874	c............................
875	c zphi is recomputed (pphi is not ok if mean molecular mass varies)
876	c with dphi = RT/mmean d(ln p) [evaluated at interface]
877
878	c DO i = 1, klon
879	c zphi(i,1) = pphis(i) + Rt_seri(i,1)/mmean(i,1)1000.
880	c * *( log(paprs(i,1)) - log(pplay(i,1)) )
881	c DO k = 2, klev
882	c zphi(i,k) = zphi(i,k-1)
883	c * + R500.(t_seri(i,k)/mmean(i,k)+t_seri(i,k-1)/mmean(i,k-1))
884	c * * (log(pplay(i,k-1)) - log(pplay(i,k)))
885	c ENDDO
886	c ENDDO
887	c............................
888	c=====
889
890	c calcul du geopotentiel aux niveaux intercouches
891	c ponderation des altitudes au niveau des couches en dp/p
892
893	DO k=1,klev
894	DO i=1,klon
895	zzlay(i,k)=zphi(i,k)/RG ! [m]
896	ENDDO
897	ENDDO
898	DO i=1,klon
899	zzlev(i,1)=pphis(i)/RG ! [m]
900	ENDDO
901	DO k=2,klev
902	DO i=1,klon
903	z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
904	z2=(paprs(i,k) +pplay(i,k))/(paprs(i,k) -pplay(i,k))
905	zzlev(i,k)=(z1zzlay(i,k-1)+z2zzlay(i,k))/(z1+z2)
906	ENDDO
907	ENDDO
908	DO i=1,klon
909	zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
910	ENDDO
911
912	c====================================================================
913	c
914	c Verifier les temperatures
915	c
916	CALL hgardfou(t_seri,ftsol,'debutphy')
917	c====================================================================
918	c
919	c Incrementer le compteur de la physique
920	c
921	itap = itap + 1
922
923	c====================================================================
924	c Orbite et eclairement
925	c====================================================================
926
927	c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0.
928	c donc pas de variations de Ls, ni de dist.
929	c La seule chose qui compte, c'est la rotation de la planete devant
930	c le Soleil...
931
932	zlongi = 0.0
933	dist = 0.72 ! en UA
934
935	c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite
936	c a sa valeur, et prendre en compte leur evolution,
937	c il faudra refaire un orbite.F...
938	c CALL orbite(zlongi,dist)
939
940	IF (cycle_diurne) THEN
941	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
942	CALL zenang(zlongi,gmtime,zdtime,latitude_deg,longitude_deg,
943	& rmu0,fract)
944	ELSE
945	call mucorr(klon,zlongi,latitude_deg,rmu0,fract)
946	ENDIF
947
948	c====================================================================
949	c Calcul des tendances traceurs
950	c====================================================================
951
952	if (iflag_trac.eq.1) then
953
954	if (tr_scheme.eq.1) then
955	! Case 1: pseudo-chemistry with relaxation toward fixed profile
956
957	call phytrac_relax (debut,lafin,nqmax,
958	I nlon,nlev,dtime,pplay,
959	O tr_seri)
960
961	elseif (tr_scheme.eq.2) then
962	! Case 2: surface emission
963	! For the moment, inspired from Mars version
964	! However, the variable 'source' could be used in physiq
965	! so the call to phytrac_emiss could be to initialise it.
966
967	call phytrac_emiss ( (rjourvrai+gmtime)*RDAY,
968	I debut,lafin,nqmax,
969	I nlon,nlev,dtime,paprs,
970	I latitude_deg,longitude_deg,
971	O tr_seri)
972
973	elseif (tr_scheme.eq.3) then ! identical to ok_chem.or.ok_cloud
974	! Case 3: Full chemistry and/or clouds
975
976	call phytrac_chimie(
977	I debut,
978	I gmtime,
979	I nqmax,
980	I klon,
981	I latitude_deg,
982	I longitude_deg,
983	I nlev,
984	I dtime,
985	I t_seri,
986	I pplay,
987	O tr_seri)
988
989	c CALL WriteField_phy('Pression',pplay,nlev)
990	c CALL WriteField_phy('PressionBnd',paprs,nlev+1)
991	c CALL WriteField_phy('Temp',t_seri,nlev)
992	c IF (ok_cloud) THEN
993	c CALL WriteField_phy('NBRTOT',NBRTOT,nlev)
994	c ENDIF
995	c CALL WriteField_phy('SAl',tr_seri(:,:,i_h2so4liq),nlev)
996	c CALL WriteField_phy('SAg',tr_seri(:,:,i_h2so4),nlev)
997
998	C === SEDIMENTATION ===
999
1000	if (ok_sedim) then
1001
1002	c !! Depend on cl_scheme !!
1003
1004	if (cl_scheme.eq.1) then
1005	c ================
1006	#ifndef MESOSCALE
1007	CALL new_cloud_sedim(
1008	I klon,
1009	I nlev,
1010	I dtime,
1011	I pplay,
1012	I paprs,
1013	I t_seri,
1014	I tr_seri,
1015	O d_tr_sed(:,:,1:2),
1016	O d_tr_ssed,
1017	I nqmax,
1018	O Fsedim)
1019
1020	DO k = 1, klev
1021	DO i = 1, klon
1022
1023	c--------------------
1024	c Ce test est necessaire pour eviter Xliq=NaN
1025	! IF (ieee_is_nan(d_tr_sed(i,k,1)).OR.
1026	! & ieee_is_nan(d_tr_sed(i,k,2))) THEN
1027	IF ((d_tr_sed(i,k,1).ne.d_tr_sed(i,k,1)).OR.
1028	& (d_tr_sed(i,k,2).ne.d_tr_sed(i,k,2))) THEN
1029	PRINT*,'sedim NaN PROBLEM'
1030	PRINT*,'d_tr_sed Nan?',d_tr_sed(i,k,:),'Temp',t_seri(i,k)
1031	PRINT*,'lat-lon',i,'level',k,'dtime',dtime
1032	PRINT*,'F_sed',Fsedim(i,k)
1033	PRINT*,'==============================================='
1034	d_tr_sed(i,k,:)=0.
1035	ENDIF
1036	c--------------------
1037
1038	tr_seri(i,k,i_h2so4liq) = tr_seri(i,k,i_h2so4liq)+
1039	& d_tr_sed(i,k,1)
1040	tr_seri(i,k,i_h2oliq) = tr_seri(i,k,i_h2oliq)+
1041	& d_tr_sed(i,k,2)
1042	d_tr_sed(i,k,:) = d_tr_sed(i,k,:) / dtime
1043	Fsedim(i,k) = Fsedim(i,k) / dtime
1044
1045	ENDDO
1046	ENDDO
1047
1048	Fsedim(:,klev+1) = 0.
1049
1050	elseif (cl_scheme.eq.2) then
1051	c ====================
1052
1053	d_tr_sed(:,:,:) = 0.D0
1054
1055	DO i=1, klon
1056
1057	c Mode 1
1058	m0_mode1(:,1)=tr_seri(i,:,i_m0_mode1drop)
1059	m0_mode1(:,2)=tr_seri(i,:,i_m0_mode1ccn)
1060	m3_mode1(:,1)=tr_seri(i,:,i_m3_mode1sa)
1061	m3_mode1(:,2)=tr_seri(i,:,i_m3_mode1w)
1062	m3_mode1(:,3)=tr_seri(i,:,i_m3_mode1ccn)
1063
1064	call drop_sedimentation(dtime,klev,m0_mode1,m3_mode1,
1065	& paprs(i,:),zzlev(i,:),zzlay(i,:),t_seri(i,:),
1066	& 1,
1067	& d_ccn_sed(:,1),d_drop_sed,d_ccn_sed(:,2),d_liq_sed)
1068
1069	d_tr_sed(i,:,i_m0_mode1drop)= d_tr_sed(i,:,i_m0_mode1drop)
1070	& + d_drop_sed
1071	d_tr_sed(i,:,i_m0_mode1ccn) = d_tr_sed(i,:,i_m0_mode1ccn)
1072	& + d_ccn_sed(:,1)
1073	d_tr_sed(i,:,i_m3_mode1ccn) = d_tr_sed(i,:,i_m3_mode1ccn)
1074	& + d_ccn_sed(:,2)
1075	d_tr_sed(i,:,i_m3_mode1sa) = d_tr_sed(i,:,i_m3_mode1sa)
1076	& + d_liq_sed(:,1)
1077	d_tr_sed(i,:,i_m3_mode1w) = d_tr_sed(i,:,i_m3_mode1w)
1078	& + d_liq_sed(:,2)
1079
1080	c Mode 2
1081	m0_mode2(:,1)=tr_seri(i,:,i_m0_mode2drop)
1082	m0_mode2(:,2)=tr_seri(i,:,i_m0_mode2ccn)
1083	m3_mode2(:,1)=tr_seri(i,:,i_m3_mode2sa)
1084	m3_mode2(:,2)=tr_seri(i,:,i_m3_mode2w)
1085	m3_mode2(:,3)=tr_seri(i,:,i_m3_mode2ccn)
1086
1087	call drop_sedimentation(dtime,klev,m0_mode2,m3_mode2,
1088	& paprs(i,:),zzlev(i,:),zzlay(i,:),t_seri(i,:),
1089	& 2,
1090	& d_ccn_sed(:,1),d_drop_sed,d_ccn_sed(:,2),d_liq_sed)
1091
1092	d_tr_sed(i,:,i_m0_mode2drop)= d_tr_sed(i,:,i_m0_mode2drop)
1093	& + d_drop_sed
1094	d_tr_sed(i,:,i_m0_mode2ccn) = d_tr_sed(i,:,i_m0_mode2ccn)
1095	& + d_ccn_sed(:,1)
1096	d_tr_sed(i,:,i_m3_mode2ccn) = d_tr_sed(i,:,i_m3_mode2ccn)
1097	& + d_ccn_sed(:,2)
1098	d_tr_sed(i,:,i_m3_mode2sa) = d_tr_sed(i,:,i_m3_mode2sa)
1099	& + d_liq_sed(:,1)
1100	d_tr_sed(i,:,i_m3_mode2w) = d_tr_sed(i,:,i_m3_mode2w)
1101	& + d_liq_sed(:,2)
1102
1103	c Aer
1104	call aer_sedimentation(dtime,klev,tr_seri(i,:,i_m0_aer),
1105	& tr_seri(i,:,i_m3_aer),paprs(i,:),
1106	& zzlev(i,:),zzlay(i,:),t_seri(i,:),
1107	& d_tr_sed(i,:,i_m0_aer),d_tr_sed(i,:,i_m3_aer),aer_flux)
1108
1109	END DO
1110
1111	DO iq=nqmax-nmicro+1,nqmax
1112	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_sed(:,:,iq)
1113	d_tr_sed(:,:,iq) = d_tr_sed(:,:,iq) / dtime
1114	END DO
1115	#endif
1116	endif
1117	c ====================
1118
1119	C === FIN SEDIMENTATION ===
1120
1121	endif ! ok_sedim
1122
1123	endif ! tr_scheme
1124	endif ! iflag_trac
1125
1126	c
1127	c====================================================================
1128	c Appeler la diffusion verticale (programme de couche limite)
1129	c====================================================================
1130
1131	c-------------------------------
1132	c VENUS TEST: on ne tient pas compte des calculs de clmain mais on force
1133	c l'equilibre radiatif du sol
1134	if (.not. ok_clmain) then
1135	if (debut) then
1136	print*,"ATTENTION, CLMAIN SHUNTEE..."
1137	endif
1138
1139	DO i = 1, klon
1140	sens(i) = 0.0e0 ! flux de chaleur sensible au sol
1141	fder(i) = 0.0e0
1142	dlw(i) = 0.0e0
1143	ENDDO
1144
1145	c Incrementer la temperature du sol
1146	c
1147	DO i = 1, klon
1148	d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
1149	ftsol(i) = ftsol(i) + d_ts(i)
1150	do j=1,nsoilmx
1151	ftsoil(i,j)=ftsol(i)
1152	enddo
1153	ENDDO
1154
1155	c-------------------------------
1156	else
1157	c-------------------------------
1158
1159	fder = dlw
1160
1161	! ADAPTATION GCM POUR CP(T)
1162
1163	if (physideal) then
1164	CALL clmain_ideal(dtime,itap,
1165	e t_seri,u_seri,v_seri,
1166	e rmu0,
1167	e ftsol,
1168	$ ftsoil,
1169	$ paprs,pplay,ppk,radsol,falbe,
1170	e solsw, sollw, sollwdown, fder,
1171	e longitude_deg, latitude_deg, dx, dy,
1172	e debut, lafin,
1173	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
1174	s fluxt,fluxu,fluxv,cdragh,cdragm,
1175	s dsens,
1176	s ycoefh,yu1,yv1)
1177	else
1178	CALL clmain(dtime,itap,
1179	e t_seri,u_seri,v_seri,
1180	e rmu0,
1181	e ftsol,
1182	$ ftsoil,
1183	$ paprs,pplay,ppk,radsol,falbe,
1184	e solsw, sollw, sollwdown, fder,
1185	e longitude_deg, latitude_deg, dx, dy,
1186	e debut, lafin,
1187	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
1188	s fluxt,fluxu,fluxv,cdragh,cdragm,
1189	s dsens,
1190	s ycoefh,yu1,yv1)
1191	endif
1192
1193	CXXX Incrementation des flux
1194	DO i = 1, klon
1195	sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
1196	fder(i) = dlw(i) + dsens(i)
1197	ENDDO
1198	CXXX
1199	IF (.not. turb_resolved) then !True only for LES
1200	DO k = 1, klev
1201	DO i = 1, klon
1202	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1203	d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s
1204	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1205	d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s
1206	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1207	d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s
1208	ENDDO
1209	ENDDO
1210	ENDIF
1211	C TRACEURS
1212
1213	if (iflag_trac.eq.1) then
1214	DO k = 1, klev
1215	DO i = 1, klon
1216	delp(i,k) = paprs(i,k)-paprs(i,k+1)
1217	ENDDO
1218	ENDDO
1219
1220	DO iq=1, nqmax
1221
1222	CALL cltrac(dtime,ycoefh,t_seri,
1223	s tr_seri(:,:,iq),source(:,iq),
1224	e paprs, pplay,delp,
1225	s d_tr_vdf(:,:,iq))
1226
1227	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
1228	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
1229
1230	ENDDO !nqmax
1231
1232	endif
1233
1234	IF (if_ebil.ge.2) THEN
1235	ztit='after clmain'
1236	CALL diagetpq(cell_area,ztit,ip_ebil,2,1,dtime
1237	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1238	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1239	call diagphy(cell_area,ztit,ip_ebil
1240	e , zero_v, zero_v, zero_v, zero_v, sens
1241	e , zero_v, zero_v, zero_v, ztsol
1242	e , d_h_vcol, d_qt, d_ec
1243	s , fs_bound, fq_bound )
1244	END IF
1245	C
1246	c
1247	c Incrementer la temperature du sol
1248	c
1249	DO i = 1, klon
1250	ftsol(i) = ftsol(i) + d_ts(i)
1251	ENDDO
1252
1253	c Calculer la derive du flux infrarouge
1254	c
1255	DO i = 1, klon
1256	dlw(i) = - 4.0RSIGMAftsol(i)**3
1257	ENDDO
1258
1259	c-------------------------------
1260	endif ! fin du VENUS TEST
1261
1262	! tests: output tendencies
1263	! call writefield_phy('physiq_d_t_vdf',d_t_vdf,klev)
1264	! call writefield_phy('physiq_d_u_vdf',d_u_vdf,klev)
1265	! call writefield_phy('physiq_d_v_vdf',d_v_vdf,klev)
1266	! call writefield_phy('physiq_d_ts',d_ts,1)
1267
1268	c
1269	c Appeler l'ajustement sec
1270	c
1271	c===================================================================
1272	c Convection seche
1273	c===================================================================
1274	c
1275	d_t_ajs(:,:)=0.
1276	d_u_ajs(:,:)=0.
1277	d_v_ajs(:,:)=0.
1278	d_tr_ajs(:,:,:)=0.
1279	c
1280	IF(prt_level>9)WRITE(lunout,*)
1281	. 'AVANT LA CONVECTION SECHE , iflag_ajs='
1282	s ,iflag_ajs
1283
1284	if(iflag_ajs.eq.0) then
1285	c Rien
1286	c ====
1287	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
1288
1289	else if(iflag_ajs.eq.1) then
1290
1291	c Ajustement sec
1292	c ==============
1293	IF(prt_level>9)WRITE(lunout,*)'ajsec'
1294
1295	! ADAPTATION GCM POUR CP(T)
1296	CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
1297	. tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
1298
1299	! ADAPTATION GCM POUR CP(T)
1300	do i=1,klon
1301	flux_ajs(i,1) = 0.0
1302	do j=2,klev
1303	flux_ajs(i,j) = flux_ajs(i,j-1)
1304	. + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime
1305	. *(paprs(i,j-1)-paprs(i,j))
1306	enddo
1307	enddo
1308
1309	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
1310	d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s
1311	u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
1312	d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s
1313	v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
1314	d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s
1315
1316	if (iflag_trac.eq.1) then
1317	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
1318	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
1319	endif
1320	endif
1321
1322	! tests: output tendencies
1323	! call writefield_phy('physiq_d_t_ajs',d_t_ajs,klev)
1324	! call writefield_phy('physiq_d_u_ajs',d_u_ajs,klev)
1325	! call writefield_phy('physiq_d_v_ajs',d_v_ajs,klev)
1326	c
1327	IF (if_ebil.ge.2) THEN
1328	ztit='after dry_adjust'
1329	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
1330	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1331	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1332	call diagphy(cell_area,ztit,ip_ebil
1333	e , zero_v, zero_v, zero_v, zero_v, sens
1334	e , zero_v, zero_v, zero_v, ztsol
1335	e , d_h_vcol, d_qt, d_ec
1336	s , fs_bound, fq_bound )
1337	END IF
1338
1339
1340	c====================================================================
1341	c RAYONNEMENT
1342	c====================================================================
1343
1344	c------------------------------------
1345	c . Compute radiative tendencies :
1346	c------------------------------------
1347	c====================================================================
1348	IF (MOD(itaprad,radpas).EQ.0) THEN
1349	c====================================================================
1350
1351	dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1352	c PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
1353
1354
1355	c------------------------------------
1356	c . Compute mean mass, cp and R :
1357	c------------------------------------
1358
1359	if(callthermos) then
1360	call concentrations2(pplay,t_seri,d_t,tr_seri, nqmax)
1361
1362	endif
1363
1364
1365	cc!!! ADD key callhedin
1366
1367	IF(callnlte.or.callthermos) THEN
1368	call compo_hedin83_mod(pplay,rmu0,
1369	& co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
1370
1371	IF(ok_chem) then
1372
1373	CC !! GG : Using only mayor species tracers abundances to compute NLTE heating/cooling
1374
1375	CC Conversion [mmr] ---> [vmr]
1376
1377	co2vmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_co2)*
1378	& mmean(1:nlon,1:nlev)/M_tr(i_co2)
1379	covmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_co)*
1380	& mmean(1:nlon,1:nlev)/M_tr(i_co)
1381	ovmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_o)*
1382	& mmean(1:nlon,1:nlev)/M_tr(i_o)
1383	n2vmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_n2)*
1384	& mmean(1:nlon,1:nlev)/M_tr(i_n2)
1385
1386	ENDIF
1387
1388	ENDIF
1389
1390	c
1391	c NLTE cooling from CO2 emission
1392	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1393
1394	IF(callnlte) THEN
1395	if(nltemodel.eq.0.or.nltemodel.eq.1) then
1396	CALL nltecool(klon, klev, nqmax, pplay*9.869e-6, t_seri,
1397	$ tr_seri, d_t_nlte)
1398	else if(nltemodel.eq.2) then
1399	CALL nlte_tcool(klon,klev,pplay*9.869e-6,
1400	$ t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
1401	$ ovmr_gcm,d_t_nlte,ierr_nlte,varerr )
1402	if(ierr_nlte.gt.0) then
1403	write(,)
1404	$ 'WARNING: nlte_tcool output with error message',
1405	$ 'ierr_nlte=',ierr_nlte,'varerr=',varerr
1406	write(,)'I will continue anyway'
1407	endif
1408
1409	endif
1410
1411	ELSE
1412
1413	d_t_nlte(:,:)=0.
1414
1415	ENDIF
1416
1417	c Find number of layers for LTE radiation calculations
1418
1419	IF(callnlte .or. callnirco2)
1420	$ CALL nlthermeq(klon, klev, paprs, pplay)
1421
1422	c
1423	c LTE radiative transfert / solar / IR matrix
1424	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1425	if (physideal) then
1426	CALL radlwsw_newtoncool
1427	e (dist, rmu0, fract, zzlev,
1428	e paprs, pplay,ftsol, t_seri)
1429	else
1430	CALL radlwsw
1431	e (dist, rmu0, fract, zzlev,
1432	e paprs, pplay,ftsol, t_seri)
1433	endif
1434
1435	c albedo variations: test for Yeon Joo Lee
1436	c increment to increase it for 20 Vd => +80%
1437	c heat(:,:)=heat(:,:)(1.+0.80((rjourvrai-356)+gmtime)/20.)
1438	c or to decrease it for 20 Vd => 1/1.8
1439	c heat(:,:)=heat(:,:)/(1.+0.80*((rjourvrai-356)+gmtime)/20.)
1440
1441	c CO2 near infrared absorption
1442	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1443
1444	d_t_nirco2(:,:)=0.
1445	if (callnirco2) then
1446	call nirco2abs (klon, klev, pplay, dist, nqmax, tr_seri,
1447	. rmu0, fract, d_t_nirco2)
1448	endif
1449
1450
1451	c Net atmospheric radiative heating rate (K.s-1)
1452	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1453
1454	IF(callnlte.or.callnirco2) THEN
1455	CALL blendrad(klon, klev, pplay,heat,
1456	& cool, d_t_nirco2,d_t_nlte, dtsw, dtlw)
1457	ELSE
1458	dtsw(:,:)=heat(:,:)
1459	dtlw(:,:)=-1*cool(:,:)
1460	ENDIF
1461
1462	DO k=1,klev
1463	DO i=1,klon
1464	d_t_rad(i,k) = dtsw(i,k) + dtlw(i,k) ! K/s
1465	ENDDO
1466	ENDDO
1467
1468
1469	cc---------------------------------------------
1470
1471	c EUV heating rate (K.s-1)
1472	c ~~~~~~~~~~~~~~~~~~~~~~~~
1473
1474	d_t_euv(:,:)=0.
1475
1476	IF (callthermos) THEN
1477
1478	c call euvheat(klon, klev,t_seri,paprs,pplay,zzlay,
1479	c $ rmu0,pdtphys,gmtime,rjourvrai, co2vmr_gcm, n2vmr_gcm,
1480	c $ covmr_gcm, ovmr_gcm,d_t_euv )
1481	call euvheat(klon, klev, nqmax, t_seri,paprs,pplay,zzlay,
1482	$ rmu0,pdtphys,gmtime,rjourvrai,
1483	$ tr_seri, d_tr, d_t_euv )
1484
1485	DO k=1,klev
1486	DO ig=1,klon
1487	d_t_rad(ig,k)=d_t_rad(ig,k)+d_t_euv(ig,k)
1488
1489	ENDDO
1490	ENDDO
1491
1492	ENDIF ! callthermos
1493
1494	c====================================================================
1495	itaprad = 0
1496	ENDIF ! radpas
1497	c====================================================================
1498	c
1499	c Ajouter la tendance des rayonnements (tous les pas)
1500	c
1501	DO k = 1, klev
1502	DO i = 1, klon
1503	t_seri(i,k) = t_seri(i,k) + d_t_rad(i,k) * dtime
1504	ENDDO
1505	ENDDO
1506
1507	! CONDUCTION and MOLECULAR VISCOSITY
1508	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1509
1510	d_t_conduc(:,:)=0.
1511	d_u_molvis(:,:)=0.
1512	d_v_molvis(:,:)=0.
1513
1514	IF (callthermos) THEN
1515
1516	tsurf(:)=t_seri(:,1)
1517	call conduction(klon, klev,pdtphys,
1518	$ pplay,paprs,t_seri,
1519	$ tsurf,zzlev,zzlay,d_t_conduc)
1520
1521	call molvis(klon, klev,pdtphys,
1522	$ pplay,paprs,t_seri,
1523	$ u,tsurf,zzlev,zzlay,d_u_molvis)
1524
1525	call molvis(klon, klev, pdtphys,
1526	$ pplay,paprs,t_seri,
1527	$ v,tsurf,zzlev,zzlay,d_v_molvis)
1528
1529	DO k=1,klev
1530	DO ig=1,klon
1531	t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K]
1532	u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! m/s
1533	v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! m/s
1534	ENDDO
1535	ENDDO
1536	ENDIF
1537
1538
1539	! -- MOLECULAR DIFFUSION ---
1540
1541	d_q_moldif(:,:,:)=0
1542
1543	IF (callthermos .and. ok_chem) THEN
1544
1545	call moldiff_red(klon, klev, nqmax,
1546	& pplay,paprs,t_seri, tr_seri, pdtphys,
1547	& zzlay,d_t_euv,d_t_conduc,d_q_moldif)
1548
1549
1550	! --- update tendencies tracers ---
1551
1552	DO iq = 1, nqmax
1553	DO k=1,klev
1554	DO ig=1,klon
1555	tr_seri(ig,k,iq)= tr_seri(ig,k,iq)+
1556	& d_q_moldif(ig,k,iq)*dtime ! [Kg/kg]?
1557	ENDDO
1558	ENDDO
1559	ENDDO
1560
1561
1562	ENDIF ! callthermos & ok_chem
1563
1564	c====================================================================
1565	! tests: output tendencies
1566	! call writefield_phy('physiq_dtrad',dtrad,klev)
1567
1568	IF (if_ebil.ge.2) THEN
1569	ztit='after rad'
1570	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
1571	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1572	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1573	call diagphy(cell_area,ztit,ip_ebil
1574	e , topsw, toplw, solsw, sollw, zero_v
1575	e , zero_v, zero_v, zero_v, ztsol
1576	e , d_h_vcol, d_qt, d_ec
1577	s , fs_bound, fq_bound )
1578	END IF
1579	c
1580
1581	c====================================================================
1582	c Calcul des gravity waves FLOTT
1583	c====================================================================
1584	c
1585	c if (ok_orodr.or.ok_gw_nonoro) then
1586
1587	c CALCUL DE N2
1588	c UTILISE LA RELATION ENTRE N2 ET STABILITE
1589	c N2 = RG/T (dT/dz+RG/cp(T))
1590	c ET DONC EN N'UTILISE QUE LA TEMPERATURE, PAS teta.
1591
1592	do i=1,klon
1593	do k=2,klev
1594	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1595	enddo
1596	enddo
1597	do i=1,klon
1598	do k=2,klev
1599	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1600	zdtlev(i,k) = t_seri(i,k)-t_seri(i,k-1)
1601	zdzlev(i,k) = (zphi(i,k)-zphi(i,k-1))/RG
1602	zn2(i,k) = RG/ztlev(i,k) * ( zdtlev(i,k)/zdzlev(i,k)
1603	. + RG/cpdet(ztlev(i,k)) )
1604	zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
1605	enddo
1606	zn2(i,1) = 1.e-12 ! securite
1607	enddo
1608
1609	c endif
1610
1611	c ----------------------------ORODRAG
1612	IF (ok_orodr) THEN
1613	c
1614	c selection des points pour lesquels le shema est actif:
1615	igwd=0
1616	DO i=1,klon
1617	itest(i)=0
1618	c IF ((zstd(i).gt.10.0)) THEN
1619	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1620	itest(i)=1
1621	igwd=igwd+1
1622	idx(igwd)=i
1623	ENDIF
1624	ENDDO
1625	c igwdim=MAX(1,igwd)
1626	c
1627	c A ADAPTER POUR VENUS!!! [ TN: c'est fait ! ]
1628	CALL drag_noro(klon,klev,dtime,paprs,pplay,pphi,zn2,
1629	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
1630	e igwd,idx,itest,
1631	e t_seri, u_seri, v_seri,
1632	s zulow, zvlow, zustrdr, zvstrdr,
1633	s d_t_oro, d_u_oro, d_v_oro,
1634	s zublstrdr,zvblstrdr,znlow,zeff,zbl,
1635	s ztau,tau0,knu2,kbreak)
1636
1637	c print*,"d_u_oro=",d_u_oro(klon/2,:)
1638	c ajout des tendances
1639	t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
1640	d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s
1641	u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
1642	d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s
1643	v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
1644	d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s
1645	c
1646	ELSE
1647	d_t_oro = 0.
1648	d_u_oro = 0.
1649	d_v_oro = 0.
1650	zustrdr = 0.
1651	zvstrdr = 0.
1652	zublstrdr = 0.
1653	zvblstrdr = 0.
1654	znlow = 0.
1655	zeff = 0.
1656	zbl = 0
1657	knu2 = 0
1658	kbreak = 0
1659	ztau = 0
1660	tau0 = 0.
1661	c
1662	ENDIF ! fin de test sur ok_orodr
1663	c
1664	! tests: output tendencies
1665	! call writefield_phy('physiq_d_t_oro',d_t_oro,klev)
1666	! call writefield_phy('physiq_d_u_oro',d_u_oro,klev)
1667	! call writefield_phy('physiq_d_v_oro',d_v_oro,klev)
1668
1669	c ----------------------------OROLIFT
1670	IF (ok_orolf) THEN
1671	print*,"ok_orolf NOT IMPLEMENTED !"
1672	stop
1673	c
1674	c selection des points pour lesquels le shema est actif:
1675	igwd=0
1676	DO i=1,klon
1677	itest(i)=0
1678	IF ((zpic(i)-zmea(i)).GT.100.) THEN
1679	itest(i)=1
1680	igwd=igwd+1
1681	idx(igwd)=i
1682	ENDIF
1683	ENDDO
1684	c igwdim=MAX(1,igwd)
1685	c
1686	c A ADAPTER POUR VENUS!!!
1687	c CALL lift_noro(klon,klev,dtime,paprs,pplay,
1688	c e latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval,
1689	c e igwd,idx,itest,
1690	c e t_seri, u_seri, v_seri,
1691	c s zulow, zvlow, zustrli, zvstrli,
1692	c s d_t_lif, d_u_lif, d_v_lif )
1693
1694	c
1695	c ajout des tendances
1696	t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
1697	d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s
1698	u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
1699	d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s
1700	v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
1701	d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s
1702	c
1703	ELSE
1704	d_t_lif = 0.
1705	d_u_lif = 0.
1706	d_v_lif = 0.
1707	zustrli = 0.
1708	zvstrli = 0.
1709	c
1710	ENDIF ! fin de test sur ok_orolf
1711
1712	c ---------------------------- NON-ORO GRAVITY WAVES
1713	IF(ok_gw_nonoro) then
1714
1715	call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
1716	e t_seri, u_seri, v_seri, paprs(klon/2+1,:),
1717	o zustrhi,zvstrhi,
1718	o d_t_hin, d_u_hin, d_v_hin)
1719
1720	c ajout des tendances
1721
1722	t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
1723	d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s
1724	u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
1725	d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s
1726	v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
1727	d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s
1728
1729	ELSE
1730	d_t_hin = 0.
1731	d_u_hin = 0.
1732	d_v_hin = 0.
1733	zustrhi = 0.
1734	zvstrhi = 0.
1735
1736	ENDIF ! fin de test sur ok_gw_nonoro
1737
1738	! tests: output tendencies
1739	! call writefield_phy('physiq_d_t_hin',d_t_hin,klev)
1740	! call writefield_phy('physiq_d_u_hin',d_u_hin,klev)
1741	! call writefield_phy('physiq_d_v_hin',d_v_hin,klev)
1742
1743	c====================================================================
1744	c Transport de ballons
1745	c====================================================================
1746	if (ballons.eq.1) then
1747	CALL ballon(30,pdtphys,rjourvrai,gmtime*RDAY,
1748	& latitude_deg,longitude_deg,
1749	c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM)
1750	C t,pplay,u,v,zphi) ! alt above planet average radius
1751	endif !ballons
1752
1753	c====================================================================
1754	c Bilan de mmt angulaire
1755	c====================================================================
1756	if (bilansmc.eq.1) then
1757	CMODDEB FLOTT
1758	C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
1759	C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
1760
1761	DO i = 1, klon
1762	zustrph(i)=0.
1763	zvstrph(i)=0.
1764	zustrcl(i)=0.
1765	zvstrcl(i)=0.
1766	ENDDO
1767	DO k = 1, klev
1768	DO i = 1, klon
1769	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
1770	c (paprs(i,k)-paprs(i,k+1))/rg
1771	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
1772	c (paprs(i,k)-paprs(i,k+1))/rg
1773	zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
1774	c (paprs(i,k)-paprs(i,k+1))/rg
1775	zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
1776	c (paprs(i,k)-paprs(i,k+1))/rg
1777	ENDDO
1778	ENDDO
1779
1780	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime*RDAY,
1781	C ra,rg,romega,
1782	C latitude_deg,longitude_deg,pphis,
1783	C zustrdr,zustrli,zustrcl,
1784	C zvstrdr,zvstrli,zvstrcl,
1785	C paprs,u,v)
1786
1787	CCMODFIN FLOTT
1788	endif !bilansmc
1789
1790	c====================================================================
1791	c====================================================================
1792	c Calculer le transport de l'eau et de l'energie (diagnostique)
1793	c
1794	c A REVOIR POUR VENUS...
1795	c
1796	c CALL transp (paprs,ftsol,
1797	c e t_seri, q_seri, u_seri, v_seri, zphi,
1798	c s ve, vq, ue, uq)
1799	c
1800	c====================================================================
1801	c+jld ec_conser
1802	DO k = 1, klev
1803	DO i = 1, klon
1804	d_t_ec(i,k)=0.5/cpdet(t_seri(i,k))
1805	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
1806	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
1807	d_t_ec(i,k) = d_t_ec(i,k)/dtime
1808	END DO
1809	END DO
1810	do i=1,klon
1811	flux_ec(i,1) = 0.0
1812	do j=2,klev
1813	flux_ec(i,j) = flux_ec(i,j-1)
1814	. +cpdet(t_seri(i,j-1))/RGd_t_ec(i,j-1)(paprs(i,j-1)-paprs(i,j))
1815	enddo
1816	enddo
1817
1818	c-jld ec_conser
1819	c====================================================================
1820	IF (if_ebil.ge.1) THEN
1821	ztit='after physic'
1822	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
1823	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1824	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1825	C Comme les tendances de la physique sont ajoute dans la dynamique,
1826	C on devrait avoir que la variation d'entalpie par la dynamique
1827	C est egale a la variation de la physique au pas de temps precedent.
1828	C Donc la somme de ces 2 variations devrait etre nulle.
1829	call diagphy(cell_area,ztit,ip_ebil
1830	e , topsw, toplw, solsw, sollw, sens
1831	e , zero_v, zero_v, zero_v, ztsol
1832	e , d_h_vcol, d_qt, d_ec
1833	s , fs_bound, fq_bound )
1834	C
1835	d_h_vcol_phy=d_h_vcol
1836	C
1837	END IF
1838	C
1839	c=======================================================================
1840	c SORTIES
1841	c=======================================================================
1842
1843	c Convertir les incrementations en tendances
1844	c
1845	DO k = 1, klev
1846	DO i = 1, klon
1847	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
1848	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
1849	d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
1850	ENDDO
1851	ENDDO
1852	c
1853	DO iq = 1, nqmax
1854	DO k = 1, klev
1855	DO i = 1, klon
1856	d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
1857	ENDDO
1858	ENDDO
1859	ENDDO
1860
1861	c------------------------
1862	c Calcul moment cinetique
1863	c------------------------
1864	c TEST VENUS...
1865	c mangtot = 0.0
1866	c DO k = 1, klev
1867	c DO i = 1, klon
1868	c mang(i,k) = RA*cos(latitude(i))
1869	c . (u_seri(i,k)+RAcos(latitude(i))*ROMEGA)
1870	c . cell_area(i)(paprs(i,k)-paprs(i,k+1))/RG
1871	c mangtot=mangtot+mang(i,k)
1872	c ENDDO
1873	c ENDDO
1874	c print*,"Moment cinetique total = ",mangtot
1875	c
1876	c------------------------
1877	c
1878	c Sauvegarder les valeurs de t et u a la fin de la physique:
1879	c
1880	DO k = 1, klev
1881	DO i = 1, klon
1882	u_ancien(i,k) = u_seri(i,k)
1883	t_ancien(i,k) = t_seri(i,k)
1884	ENDDO
1885	ENDDO
1886	c
1887	c=============================================================
1888	c Ecriture des sorties
1889	c=============================================================
1890	#ifndef MESOSCALE
1891	#ifdef CPP_IOIPSL
1892
1893	#ifdef histhf
1894	#include "write_histhf.h"
1895	#endif
1896
1897	#ifdef histday
1898	#include "write_histday.h"
1899	#endif
1900
1901	#ifdef histmth
1902	#include "write_histmth.h"
1903	#endif
1904
1905	#ifdef histins
1906	#include "write_histins.h"
1907	#endif
1908
1909	#endif
1910
1911	! XIOS outputs
1912	! This can be done ANYWHERE in the physics routines !
1913
1914	#ifdef CPP_XIOS
1915	! Send fields to XIOS: (NB these fields must also be defined as
1916	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
1917
1918	! 2D fields
1919
1920	CALL send_xios_field("phis",pphis)
1921	cell_area_out(:)=cell_area(:)
1922	if (is_north_pole_phy) cell_area_out(1)=cell_area(1)/nbp_lon
1923	if (is_south_pole_phy) cell_area_out(klon)=cell_area(klon)/nbp_lon
1924	CALL send_xios_field("aire",cell_area_out)
1925	CALL send_xios_field("tsol",ftsol)
1926	CALL send_xios_field("psol",paprs(:,1))
1927	CALL send_xios_field("cdragh",cdragh)
1928	CALL send_xios_field("cdragm",cdragm)
1929
1930	CALL send_xios_field("tops",topsw)
1931	CALL send_xios_field("topl",toplw)
1932	CALL send_xios_field("sols",solsw)
1933	CALL send_xios_field("soll",sollw)
1934
1935	! 3D fields
1936
1937	CALL send_xios_field("temp",t_seri)
1938	CALL send_xios_field("pres",pplay)
1939	CALL send_xios_field("geop",zphi)
1940	CALL send_xios_field("vitu",u_seri)
1941	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
1942	CALL send_xios_field("vitv",-1.*v_seri)
1943	CALL send_xios_field("vitw",omega)
1944	CALL send_xios_field("Kz",ycoefh)
1945	CALL send_xios_field("mmean",mmean)
1946	CALL send_xios_field("rho",rho)
1947	CALL send_xios_field("BV2",zn2)
1948
1949	CALL send_xios_field("dudyn",d_u_dyn)
1950	CALL send_xios_field("duvdf",d_u_vdf)
1951	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
1952	CALL send_xios_field("dvvdf",-1.*d_v_vdf)
1953	CALL send_xios_field("duajs",d_u_ajs)
1954	CALL send_xios_field("dugwo",d_u_oro)
1955	CALL send_xios_field("dugwno",d_u_hin)
1956	CALL send_xios_field("dumolvis",d_u_molvis)
1957	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
1958	CALL send_xios_field("dvmolvis",-1.*d_v_molvis)
1959	CALL send_xios_field("dtdyn",d_t_dyn)
1960	CALL send_xios_field("dtphy",d_t)
1961	CALL send_xios_field("dtvdf",d_t_vdf)
1962	CALL send_xios_field("dtajs",d_t_ajs)
1963	CALL send_xios_field("dtswr",dtsw)
1964	CALL send_xios_field("dtswrNLTE",d_t_nirco2)
1965	CALL send_xios_field("dtswrLTE",heat)
1966	CALL send_xios_field("dtlwr",dtlw)
1967	CALL send_xios_field("dtlwrNLTE",d_t_nlte)
1968	CALL send_xios_field("dtlwrLTE",-1.*cool)
1969	CALL send_xios_field("dteuv",d_t_euv)
1970	CALL send_xios_field("dtcond",d_t_conduc)
1971	CALL send_xios_field("dtec",d_t_ec)
1972
1973	CALL send_xios_field("SWnet",swnet(:,1:klev))
1974	CALL send_xios_field("LWnet",lwnet(:,1:klev))
1975	CALL send_xios_field("fluxvdf",fluxt)
1976	CALL send_xios_field("fluxdyn",flux_dyn)
1977	CALL send_xios_field("fluxajs",flux_ajs)
1978	CALL send_xios_field("fluxec",flux_ec)
1979
1980	c When using tracers
1981	IF (iflag_trac.eq.1) THEN
1982	c photochemical compounds !!!outputs in [vmr]
1983	DO iq=1,nqmax-nmicro
1984	CALL send_xios_field(tname(iq),qx(:,:,iq)*mmean(:,:)/M_tr(iq))
1985	ENDDO
1986
1987	IF ((tr_scheme.eq.3).and.(cl_scheme.eq.1)) THEN
1988	c liquids !!!outputs in [vmr]
1989	CALL send_xios_field(tname(i_h2oliq),
1990	. qx(:,:,i_h2oliq)*mmean(:,:)/M_tr(i_h2oliq))
1991	CALL send_xios_field(tname(i_h2so4liq),
1992	. qx(:,:,i_h2so4liq)*mmean(:,:)/M_tr(i_h2so4liq))
1993	if (ok_sedim) CALL send_xios_field("Fsedim",Fsedim(:,1:klev))
1994	ENDIF
1995	ENDIF
1996
1997	IF (callthermos .and. ok_chem) THEN
1998	CALL send_xios_field("d_qmoldifCO2",d_q_moldif(:,:,i_co2))
1999	CALL send_xios_field("d_qmoldifO3p",d_q_moldif(:,:,i_o))
2000	CALL send_xios_field("d_qmoldifN2",d_q_moldif(:,:,i_n2))
2001	ENDIF
2002
2003	if (lafin.and.is_omp_master) then
2004	write(,) "physiq: call xios_context_finalize"
2005	call xios_context_finalize
2006	endif
2007
2008	#endif
2009	#else
2010	! Outputs MESOSCALE
2011	CALL allocate_comm_wrf(klon,klev)
2012	comm_HR_SW(1:klon,1:klev) = dtsw(1:klon,1:klev)
2013	comm_HR_LW(1:klon,1:klev) = dtlw(1:klon,1:klev)
2014	comm_DT_RAD(1:klon,1:klev) = d_t_rad(1:klon,1:klev)
2015	IF (turb_resolved) THEN
2016	open(17,file='hrdyn.txt',form='formatted',status='old')
2017	rewind(17)
2018	DO k=1,klev
2019	read(17,*) dt_dyn(k)
2020	ENDDO
2021	close(17)
2022
2023	do i=1,klon
2024	d_t(i,:)=d_t(i,:)+dt_dyn(:)
2025	comm_HR_DYN(i,:) = dt_dyn(:)
2026	enddo
2027	ELSE
2028	comm_HR_DYN(1:klon,1:klev) = d_t_dyn(1:klon,1:klev)
2029	comm_DT_VDF(1:klon,1:klev) = d_t_vdf(1:klon,1:klev)
2030	comm_DT_AJS(1:klon,1:klev) = d_t_ajs(1:klon,1:klev)
2031	ENDIF
2032	comm_DT(1:klon,1:klev)=d_t(1:klon,1:klev)
2033	#endif
2034
2035
2036	c====================================================================
2037	c Si c'est la fin, il faut conserver l'etat de redemarrage
2038	c====================================================================
2039	c
2040	IF (lafin) THEN
2041	itau_phy = itau_phy + itap
2042	CALL phyredem ("restartphy.nc")
2043
2044	c--------------FLOTT
2045	CMODEB LOTT
2046	C FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
2047	C DU BILAN DE MOMENT ANGULAIRE.
2048	if (bilansmc.eq.1) then
2049	write(,)'Fermeture de aaam_bud.out (FL Vous parle)'
2050	close(27)
2051	close(28)
2052	endif !bilansmc
2053	CMODFIN
2054	c-------------
2055	c--------------SLEBONNOIS
2056	C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
2057	C DES BALLONS
2058	if (ballons.eq.1) then
2059	write(,)'Fermeture des ballons*.out'
2060	close(30)
2061	close(31)
2062	close(32)
2063	close(33)
2064	close(34)
2065	endif !ballons
2066	c-------------
2067	ENDIF
2068
2069	END SUBROUTINE physiq
2070
2071	END MODULE physiq_mod
2072

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

Download in other formats:

Original Format