Context Navigation

physiq.F @ 1543

Last change on this file since 1543 was 1543, checked in by emillour, 9 years ago

All models: Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation:

dyn3d:
adapted gcm.F so that all physics initializations are now done in iniphysiq.

dyn3dpar:
adapted gcm.F so that all physics initializations are now done in iniphysiq.
updated calfis_p.F to follow up with changes.
copied over updated "bands.F90" from LMDZ5.

dynphy_lonlat:
calfis_p.F90, mod_interface_dyn_phys.F90, follow up of changes in phy_common/mod_* routines

phy_common:
added "geometry_mod.F90" to store information about the grid (replaces phy*/comgeomphy.F90) and give variables friendlier names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
added "physics_distribution_mod.F90"
updated "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_mpi_data.F90", "mod_phys_lmdz_para.F90", "mod_phys_lmdz_mpi_transfert.F90", "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_omp_data.F90", "mod_phys_lmdz_omp_transfert.F90", "write_field_phy.F90" and "ioipsl_getin_p_mod.F90" to LMDZ5 versions.

phy[venus/titan/mars/std]:
removed "init_phys_lmdz.F90", "comgeomphy.F90"; adapted routines to use geometry_mod (longitude, latitude, cell_area, etc.)

File size: 53.8 KB

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

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

Download in other formats:

Original Format