Context Navigation

physiq_mod.F @ 1658

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

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

Download in other formats:

Original Format