source: trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F @ 2008

Last change on this file since 2008 was 1726, checked in by slebonnois, 8 years ago

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