Context Navigation

physiq.F @ 1524

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

All GCMS:
More updates to enforce dynamics/physics separation:

get rid of references to "temps_mod" from physics packages;
make a "time_phylmdz_mod.F90" module to store that
information and fill it via "iniphysiq".

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

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

Download in other formats:

Original Format