source: trunk/LMDZ.TITAN.old/libf/phytitan/physiq_mod.F @ 3094

Last change on this file since 3094 was 1621, checked in by emillour, 8 years ago

Further work on full dynamics/physics separation.

LMDZ.COMMON:

  • added phy_common/vertical_layers_mod.F90 to store information on vertical grid. This is where routines in the physics should get the information.
  • The contents of vertical_layers_mod intialized via dynphy_lonlat/inigeomphy_mod.F90.

LMDZ.MARS:

  • physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
  • created an "ini_tracer_mod" routine in module "tracer_mod" for a cleaner initialization of the later.
  • removed some purely dynamics-related outputs (etot0, zoom parameters, etc.) from diagfi.nc and stats.nc outputs as these informations are not available in the physics.

LMDZ.GENERIC:

  • physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
  • added nqtot to tracer_h.F90.
  • removed some purely dynamics-related outputs (etot0, zoom parameters, etc.) from diagfi.nc and stats.nc outputs as these informations are not available in the physics.

LMDZ.VENUS:

  • physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
  • added infotrac_phy.F90 to store information on tracers in the physics. Initialized via iniphysiq.
  • added cpdet_phy_mod.F90 to store t2tpot etc. functions to be used in the physics. Initialized via iniphysiq. IMPORTANT: there are some hard-coded constants! These should match what is in cpdet_mod.F90 in the dynamics.
  • got rid of references to moyzon_mod module within the physics. The required variables (tmoy, plevmoy) are passed to the physics as arguments to physiq.

LMDZ.TITAN:

  • added infotrac_phy.F90 to store information on tracers in the physics. Initialized via iniphysiq.
  • added cpdet_phy_mod.F90 to store t2tpot etc. functions to be used in the physics.
  • Extra work required to completely decouple physics and dynamics: moyzon_mod should be cleaned up and information passed from dynamics to physics as as arguments. Likewise moyzon_ch and moyzon_mu should not be queried from logic_mod (which is in the dynamics).

EM
File size: 50.1 KB
Line 
1!
2! $Id: $
3!
4      MODULE physiq_mod
5
6      IMPLICIT NONE
7
8      CONTAINS
9
10      SUBROUTINE physiq (nlon,nlev,nqmax,
11     .            debut,lafin,rjourvrai,gmtime,pdtphys,
12     .            paprs,pplay,ppk,pphi,pphis,presnivs,
13     .            u,v,t,qx,
14     .            flxmw,
15     .            d_u, d_v, d_t, d_qx, d_ps)
16
17c======================================================================
18c
19c Modifications pour la physique de Titan
20c     S. Lebonnois (LMD/CNRS) Juin 2013: Parallelisation
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-temps universel dans la journee (0 a RDAY s)
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 infotrac_phy, ONLY: iflag_trac, tname, ttext
64      use dimphy
65      USE geometry_mod, ONLY: longitude, latitude, ! in radians
66     &                        longitude_deg, latitude_deg, ! in degrees
67     &                        cell_area, dx, dy
68      use cpdet_phy_mod, only: cpdet, t2tpot
69      USE mod_phys_lmdz_para, only : is_parallel,jj_nb,
70     &                               is_north_pole_phy,
71     &                               is_south_pole_phy
72      USE phys_state_var_mod ! Variables sauvegardees de la physique
73      USE iophy
74      USE common_mod, only: rmcbar,xfbar,ncount,
75     &      flxesp_i,tau_drop,tau_aer,solesp,precip,
76     &      evapch4,occcld_m,occcld,satch4,satc2h6,satc2h2,rmcloud,
77     &      TauHID,TauHVD,TauGID,TauGVD,TauCID,TauCVD,NSPECV,NSPECI,
78     &      common_init
79
80      USE moyzon_mod
81      USE write_field_phy
82      USE time_phylmdz_mod, only: itau_phy,day_ref,annee_ref,nday
83      USE logic_mod, only: moyzon_ch,moyzon_mu
84      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
85      IMPLICIT none
86c======================================================================
87c   CLEFS CPP POUR LES IO
88c   =====================
89#define histday
90#define histmth
91#define histins
92c======================================================================
93#include "dimensions.h"
94      integer jjmp1
95      parameter (jjmp1=jjm+1-1/jjm)
96#include "dimsoil.h"
97#include "clesphys.h"
98#include "iniprint.h"
99#include "tabcontrol.h"
100#include "comorbit.h"
101#include "microtab.h"
102#include "itemps.h"
103c======================================================================
104      LOGICAL ok_journe ! sortir le fichier journalier
105      save ok_journe
106c      PARAMETER (ok_journe=.true.)
107c
108      LOGICAL ok_mensuel ! sortir le fichier mensuel
109      save ok_mensuel
110c      PARAMETER (ok_mensuel=.true.)
111c
112      LOGICAL ok_instan ! sortir le fichier instantane
113      save ok_instan
114c      PARAMETER (ok_instan=.true.)
115c
116c======================================================================
117c
118c Variables argument:
119c
120      INTEGER nlon
121      INTEGER nlev
122      INTEGER nqmax
123      REAL rjourvrai
124      REAL gmtime
125      REAL pdtphys
126      LOGICAL debut, lafin
127      REAL paprs(klon,klev+1)
128      REAL pplay(klon,klev)
129      REAL pphi(klon,klev)
130      REAL pphis(klon)
131      REAL presnivs(klev)
132
133! ADAPTATION GCM POUR CP(T)
134      REAL ppk(klon,klev)
135
136      REAL u(klon,klev)
137      REAL v(klon,klev)
138      REAL t(klon,klev)
139      REAL qx(klon,klev,nqmax)
140
141      REAL d_u_dyn(klon,klev)
142      REAL d_t_dyn(klon,klev)
143
144      REAL flxmw(klon,klev)
145
146      REAL d_u(klon,klev)
147      REAL d_v(klon,klev)
148      REAL d_t(klon,klev)
149      REAL d_qx(klon,klev,nqmax)
150      REAL d_ps(klon)
151
152c Variables propres a la physique
153c
154      REAL,save,allocatable :: rlev(:,:) ! altitude a chaque niveau (interface inferieure de la couche)
155      INTEGER,save :: itap        ! compteur pour la physique
156      REAL delp(klon,klev)        ! epaisseur d'une couche
157      REAL omega(klon,klev)
158     
159      INTEGER igwd,idx(klon),itest(klon)
160c
161c  Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro
162
163      REAL zulow(klon),zvlow(klon)
164      REAL zustrdr(klon), zvstrdr(klon)
165      REAL zustrli(klon), zvstrli(klon)
166      REAL zustrhi(klon), zvstrhi(klon)
167
168c Pour calcul GW drag oro et nonoro: CALCUL de N2:
169      real zdzlev(klon,klev)
170      real ztlev(klon,klev),zpklev(klon,klev)
171      real ztetalay(klon,klev),ztetalev(klon,klev)
172      real zdtetalev(klon,klev)
173      real zn2(klon,klev) ! BV^2 at plev
174
175c Pour les bilans de moment angulaire,
176      integer bilansmc
177c Pour le transport de ballons
178      integer ballons
179c j'ai aussi besoin
180c du stress de couche limite a la surface:
181
182      REAL zustrcl(klon),zvstrcl(klon)
183
184c et du stress total c de la physique:
185
186      REAL zustrph(klon),zvstrph(klon)
187
188c Variables locales:
189c
190      REAL cdragh(klon) ! drag coefficient pour T and Q
191      REAL cdragm(klon) ! drag coefficient pour vent
192c
193cAA  Pour  TRACEURS
194cAA
195      REAL,save,allocatable :: source(:,:)
196      integer nmicro
197      save    nmicro
198      character*8 nom
199      REAL qaer(klon,klev,nqmax)
200
201      REAL ycoefh(klon,klev)    ! coef d'echange pour phytrac
202      REAL yu1(klon)            ! vents dans la premiere couche U
203      REAL yv1(klon)            ! vents dans la premiere couche V
204
205      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
206      REAL ve(klon) ! integr. verticale du transport meri. de l'energie
207      REAL vq(klon) ! integr. verticale du transport meri. de l'eau
208      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
209      REAL uq(klon) ! integr. verticale du transport zonal de l'eau
210c
211
212c======================================================================
213c
214c Declaration des procedures appelees
215c
216      EXTERNAL ajsec     ! ajustement sec
217      EXTERNAL clmain    ! couche limite
218      EXTERNAL hgardfou  ! verifier les temperatures
219      EXTERNAL orbite    ! calculer l'orbite
220      EXTERNAL phyetat0  ! lire l'etat initial de la physique
221      EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
222      EXTERNAL radlwsw   ! rayonnements solaire et infrarouge
223      EXTERNAL suphec    ! initialiser certaines constantes
224c     EXTERNAL transp    ! transport total de l'eau et de l'energie
225      EXTERNAL abort_gcm
226      EXTERNAL printflag
227      EXTERNAL zenang
228      EXTERNAL diagetpq
229      EXTERNAL conf_phys
230      EXTERNAL diagphy
231      EXTERNAL mucorr
232      EXTERNAL phytrac
233c
234c Variables locales
235c
236CXXX PB
237      REAL fluxt(klon,klev)   ! flux turbulent de chaleur
238      REAL fluxu(klon,klev)   ! flux turbulent de vitesse u
239      REAL fluxv(klon,klev)   ! flux turbulent de vitesse v
240c
241      REAL flux_dyn(klon,klev)  ! flux de chaleur produit par la dynamique
242      REAL flux_ajs(klon,klev)  ! flux de chaleur ajustement sec
243      REAL flux_ec(klon,klev)   ! flux de chaleur Ec
244c
245      REAL    dtimerad
246      INTEGER itaprad
247      SAVE itaprad,dtimerad
248      REAL zdtime
249c
250c CHIMIE
251
252      REAL    dtimechim
253      INTEGER itapchim,appel_chim
254      SAVE itapchim,dtimechim
255
256c ORBITE
257
258      REAL dist, rmu0(klon), fract(klon), pdecli
259      REAL rmu0bar(klon), fractbar(klon)
260      REAL zday
261      REAL zls,zlsdeg,zlsm1
262      save zlsm1
263c
264      INTEGER i, k, iq, ig, j, ll, l
265c
266      REAL zphi(klon,klev)
267      REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
268c
269c Variables du changement
270c
271c ajs: ajustement sec
272c vdf: couche limite (Vertical DiFfusion)
273c mph: microphysique
274c kim: chimie
275      REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax)
276      REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
277c
278      REAL d_ts(klon)
279c
280      REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
281      REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax)
282c
283      REAL d_tr_mph(klon,klev,nqmax),d_tr_kim(klon,klev,nqmax)
284
285CMOD 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)
290C          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
294c
295c Variables liees a l'ecriture de la bande histoire physique
296c
297      INTEGER ecrit_mth
298      SAVE ecrit_mth   ! frequence d'ecriture (fichier mensuel)
299c
300      INTEGER ecrit_day
301      SAVE ecrit_day   ! frequence d'ecriture (fichier journalier)
302c
303      INTEGER ecrit_ins
304      SAVE ecrit_ins   ! frequence d'ecriture (fichier instantane)
305c
306      integer itau_w   ! pas de temps ecriture = itap + itau_phy
307
308c Variables locales pour effectuer les appels en serie
309c
310      REAL t_seri(klon,klev)
311      REAL u_seri(klon,klev), v_seri(klon,klev)
312c
313      REAL tr_seri(klon,klev,nqmax)
314      REAL d_tr(klon,klev,nqmax)
315c
316c pour ioipsl
317      INTEGER nid_day, nid_mth, nid_ins
318      SAVE nid_day, nid_mth, nid_ins
319      INTEGER nhori, nvert, idayref
320      REAL zsto, zout, zsto1, zsto2, zero
321      parameter (zero=0.0e0)
322      real zjulian
323      save zjulian
324      REAL tmpout(klon,klev)  ! pour sorties
325
326      CHARACTER*1  str1
327      CHARACTER*2  str2
328      character*20 modname
329      character*80 abort_message
330      logical ok_sync
331
332      character*30 nom_fichier
333      character*10 varname
334      character*40 vartitle
335      character*20 varunits
336C     Variables liees au bilan d'energie et d'enthalpi
337      REAL ztsol(klon)
338      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
339     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
340      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
341     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
342      REAL      d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
343      REAL      d_h_vcol_phy
344      REAL      fs_bound, fq_bound
345      SAVE      d_h_vcol_phy
346      REAL      zero_v(klon),zero_v2(klon,klev)
347      CHARACTER*15 ztit
348      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
349      SAVE      ip_ebil
350      DATA      ip_ebil/2/
351      INTEGER   if_ebil ! level for energy conserv. dignostics
352      SAVE      if_ebil
353c+jld ec_conser
354      REAL d_t_ec(klon,klev)    ! tendance du a la conversion Ec -> E thermique
355c-jld ec_conser
356
357c TEST VENUS...
358      REAL mang(klon,klev)    ! moment cinetique
359      REAL mangtot            ! moment cinetique total
360
361c Temporaire avant de trouver mieux :
362c Recuperation des TAU du TR
363      REAL t_tauhvd(klon,klev),t_khvd(klon,klev)
364      REAL t_tcld(klon,klev),t_kcld(klon,klev)
365      REAL t_kcvd(klon,klev)
366
367       REAL ch4(klon,jjm+1),dch4(jjm+1)
368       INTEGER ig0
369       integer ich4
370       common/ch4ind/ich4
371
372c     flux de chaleur latente d'evaporation CH4
373      REAL fclat(klon)
374c     reservoir de surface
375      REAL,save,allocatable :: reservoir(:)
376
377c cell_area for outputs in hist*
378      REAL cell_area_out(klon)
379
380c Declaration des constantes et des fonctions thermodynamiques
381c
382#include "YOMCST.h"
383
384c======================================================================
385c INITIALISATIONS
386c================
387
388      modname = 'physiq'
389      ok_sync=.TRUE.
390
391      bilansmc = 0
392      ballons  = 0
393! NE FONCTIONNENT PAS ENCORE EN PARALLELE !!!
394      if (is_parallel) then
395        bilansmc = 0
396        ballons  = 0
397      endif
398
399      IF (if_ebil.ge.1) THEN
400        DO i=1,klon
401          zero_v(i)=0.
402        END DO
403        DO i=1,klon
404         DO j=1,klev
405          zero_v2(i,j)=0.
406         END DO
407        END DO
408      END IF
409     
410c PREMIER APPEL SEULEMENT
411c========================
412      IF (debut) THEN
413         allocate(rlev(klon,klevp1))
414         allocate(source(klon,nqmax))
415         allocate(reservoir(klon))
416
417         CALL suphec ! initialiser constantes et parametres phys.
418
419         IF (if_ebil.ge.1) d_h_vcol_phy=0.
420c
421c appel a la lecture du physiq.def
422c
423         call conf_phys(ok_mensuel,ok_journe,
424     .                  ok_instan,
425     .                  if_ebil)
426
427         call phys_state_var_init
428         call common_init
429c
430c Initialiser les compteurs:
431c
432         itap    = 0
433         itaprad = 0
434         itapchim    = 1
435
436c init rnuabar
437         ncount(:,:) = 0
438         rmcbar  = 0.
439         xfbar   = 0.
440         
441c         
442c Lecture startphy.nc :
443c
444         CALL phyetat0 ("startphy.nc")
445
446c dtime est defini dans tabcontrol.h et lu dans startphy
447c pdtphys est calcule a partir des nouvelles conditions:
448c Reinitialisation du pas de temps physique quand changement
449         IF (ABS(dtime-pdtphys).GT.0.001) THEN
450            WRITE(lunout,*) 'Pas physique a change',dtime,
451     .                        pdtphys
452c           abort_message='Pas physique n est pas correct '
453c           call abort_gcm(modname,abort_message,1)
454c----------------
455c pour initialiser convenablement le time_counter, il faut tenir compte
456c du changement de dtime en changeant itau_phy (point de depart)
457            itau_phy = NINT(itau_phy*dtime/pdtphys)
458c----------------
459            dtime=pdtphys
460         ENDIF
461
462         radpas  = NINT( RDAY/pdtphys/nbapp_rad)
463         chimpas =   radpas*nbapp_rad/nbapp_chim
464
465         CALL printflag( ok_mensuel,ok_journe,ok_instan )
466c
467c Initialiser les pas de temps:
468c
469      dtimerad = dtime*REAL(radpas)  ! pas de temps du rayonnement (s)
470c      PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
471           
472      dtimechim = dtime*REAL(chimpas)  ! pas de temps de la chimie (s)
473c      PRINT*,'dtimechim,dtime,chimpas',dtimechim,dtime,chimpas
474
475
476c INITIALISATION ORBITE
477
478         CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
479
480c---------
481c FLOTT
482       IF (ok_orodr) THEN
483         DO i=1,klon
484         rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
485         ENDDO
486         CALL SUGWD(klon,klev,paprs,pplay)
487         DO i=1,klon
488         zuthe(i)=0.
489         zvthe(i)=0.
490         if(zstd(i).gt.10.)then
491           zuthe(i)=(1.-zgam(i))*cos(zthe(i))
492           zvthe(i)=(1.-zgam(i))*sin(zthe(i))
493         endif
494         ENDDO
495       ENDIF
496
497      if (bilansmc.eq.1) then
498C  OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
499C  DU BILAN DE MOMENT ANGULAIRE.
500      open(27,file='aaam_bud.out',form='formatted')
501      open(28,file='fields_2d.out',form='formatted')
502      write(*,*)'Ouverture de aaam_bud.out (FL Vous parle)'
503      write(*,*)'Ouverture de fields_2d.out (FL Vous parle)'
504      endif !bilansmc
505
506c--------------SLEBONNOIS
507C  OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
508C  DES BALLONS
509      if (ballons.eq.1) then
510      open(30,file='ballons-lat.out',form='formatted')
511      open(31,file='ballons-lon.out',form='formatted')
512      open(32,file='ballons-u.out',form='formatted')
513      open(33,file='ballons-v.out',form='formatted')
514      open(34,file='ballons-alt.out',form='formatted')
515      write(*,*)'Ouverture des ballons*.out'
516      endif !ballons
517c-------------
518
519c---------
520C TRACEURS
521C source dans couche limite
522         source = 0.0 ! pas de source, pour l'instant
523C
524c Si microphysique offline, pas besoin d'avoir de traceurs microphysiques
525c car on lit les profils verticaux des qaer dans une look-up table pour
526c le rayonnement.
527
528c  calcul de nmicro
529c !!!! Les traceurs microphysiques doivent etre toujours en premiers!!
530
531      nmicro = 0
532      do iq=1,nqmax
533         nom = tname(iq)
534c        print*,iq,"nom=",nom,"tname=",tname(iq)
535         print*,iq,"nom=",nom
536         if (nom(1:1).eq."q") then
537           nmicro = nmicro+1
538         endif
539      enddo
540      print*,"nmicro=",nmicro
541
542c --------------
543c Verifications:
544c --------------
545         IF ((nmicro.eq.0).and.(microfi.eq.1)) THEN
546           abort_message="MICROPHYSIQUE ONLINE, MAIS NMICRO=0..."
547           call abort_gcm(modname,abort_message,1)
548         ENDIF
549         IF (microfi.lt.1.and.clouds.eq.1) THEN
550          write(lunout,*)"microfi.lt.1.and.clouds.eq.1"
551          abort_message =
552     &    "Impossible de faire des nuages sans microphysique..."
553          call abort_gcm(modname,abort_message,1)
554         ENDIF
555         IF (nlon .NE. klon) THEN
556            WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
557     .                      klon
558            abort_message='nlon et klon ne sont pas coherents'
559            call abort_gcm(modname,abort_message,1)
560         ENDIF
561         IF (nlev .NE. klev) THEN
562            WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
563     .                       klev
564            abort_message='nlev et klev ne sont pas coherents'
565            call abort_gcm(modname,abort_message,1)
566         ENDIF
567
568         IF (((moyzon_mu).and.(microfi.ne.1)).or.
569     .       ((.not.moyzon_mu).and.(microfi.eq.1))) THEN
570           abort_message="Microphysic 2D and moyzon_mu not compatible"
571           write(lunout,*) "moyzon_mu=",moyzon_mu
572           write(lunout,*) "microfi=",microfi
573           call abort_gcm(modname,abort_message,1)
574         ENDIF
575         IF (((moyzon_ch).and.(.not.chimi)).or.
576     .       ((.not.moyzon_ch).and.(chimi))) THEN
577           abort_message="Chemistry and moyzon_ch not compatible"
578           write(lunout,*) "moyzon_ch=",moyzon_ch
579           write(lunout,*) "chimi=",chimi
580           call abort_gcm(modname,abort_message,1)
581         ENDIF
582
583         IF (dtime*REAL(radpas).GT.(RDAY*0.25).AND.cycle_diurne)
584     $    THEN
585           WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
586           WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
587           abort_message='Nbre d appels au rayonnement insuffisant'
588           call abort_gcm(modname,abort_message,1)
589         ENDIF
590c
591         WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
592     .                   iflag_ajs
593c
594         ecrit_mth = NINT(RDAY/dtime) *nday  ! tous les nday jours
595         IF (ok_mensuel) THEN
596         WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
597     .                   ecrit_mth
598         ENDIF
599
600         ecrit_day = NINT(RDAY/dtime *1.0)  ! tous les jours
601         IF (ok_journe) THEN
602         WRITE(lunout,*)'La frequence de sortie journaliere est de ',
603     .                   ecrit_day
604         ENDIF
605
606         ecrit_ins = NINT(RDAY/dtime*ecriphy)  ! Fraction de jour reglable
607         IF (ok_instan) THEN
608         WRITE(lunout,*)'La frequence de sortie instant. est de ',
609     .                   ecrit_ins
610         ENDIF
611
612c Initialisation des sorties
613c===========================
614
615#ifdef CPP_IOIPSL
616
617#ifdef histday
618#include "ini_histday.h"
619#endif
620
621#ifdef histmth
622#include "ini_histmth.h"
623#endif
624
625#ifdef histins
626#include "ini_histins.h"
627#endif
628
629#endif
630
631c
632c Initialiser les valeurs de u pour calculs tendances
633c (pour T, c'est fait dans phyetat0)
634c
635      DO k = 1, klev
636      DO i = 1, klon
637         u_ancien(i,k) = u(i,k)
638      ENDDO
639      ENDDO
640
641      ENDIF ! debut
642c====================================================================
643c======================================================================
644
645c   Creer un reservoir de surface infini
646c
647      reservoir(:) = 2.
648
649c Mettre a zero des variables de sortie (pour securite)
650c
651      DO i = 1, klon
652         d_ps(i) = 0.0
653      ENDDO
654      DO k = 1, klev
655      DO i = 1, klon
656         d_t(i,k) = 0.0
657         d_u(i,k) = 0.0
658         d_v(i,k) = 0.0
659      ENDDO
660      ENDDO
661      DO iq = 1, nqmax
662      DO k = 1, klev
663      DO i = 1, klon
664         d_qx(i,k,iq) = 0.0
665      ENDDO
666      ENDDO
667      ENDDO
668c
669c Ne pas affecter les valeurs entrees de u, v, h, et q
670c
671      DO k = 1, klev
672      DO i = 1, klon
673         t_seri(i,k)  = t(i,k)
674         u_seri(i,k)  = u(i,k)
675         v_seri(i,k)  = v(i,k)
676      ENDDO
677      ENDDO
678      DO iq = 1, nqmax
679      DO  k = 1, klev
680      DO  i = 1, klon
681         tr_seri(i,k,iq) = qx(i,k,iq)
682      ENDDO
683      ENDDO
684      ENDDO
685C
686      DO i = 1, klon
687          ztsol(i) = ftsol(i)
688      ENDDO
689C
690      IF (if_ebil.ge.1) THEN
691        ztit='after dynamic'
692        CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
693     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
694     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
695C     Comme les tendances de la physique sont ajoute dans la dynamique,
696C     on devrait avoir que la variation d'entalpie par la dynamique
697C     est egale a la variation de la physique au pas de temps precedent.
698C     Donc la somme de ces 2 variations devrait etre nulle.
699        call diagphy(cell_area,ztit,ip_ebil
700     e      , zero_v, zero_v, zero_v, zero_v, zero_v
701     e      , zero_v, zero_v, zero_v, ztsol
702     e      , d_h_vcol+d_h_vcol_phy, d_qt, 0.
703     s      , fs_bound, fq_bound )
704      END IF
705
706c====================================================================
707c Diagnostiquer la tendance dynamique
708c
709      IF (ancien_ok) THEN
710         DO k = 1, klev
711         DO i = 1, klon
712            d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
713            d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
714         ENDDO
715         ENDDO
716
717! ADAPTATION GCM POUR CP(T)
718         do i=1,klon
719          flux_dyn(i,1) = 0.0
720          do j=2,klev
721            flux_dyn(i,j) = flux_dyn(i,j-1)
722     . +cpdet(t_seri(i,j-1))/RG*d_t_dyn(i,j-1)*(paprs(i,j-1)-paprs(i,j))
723          enddo
724         enddo
725         
726      ELSE
727         DO k = 1, klev
728         DO i = 1, klon
729            d_u_dyn(i,k) = 0.0
730            d_t_dyn(i,k) = 0.0
731         ENDDO
732         ENDDO
733         ancien_ok = .TRUE.
734      ENDIF
735c====================================================================
736c
737c Calcule de vitesse verticale a partir de flux de masse verticale
738      DO k = 1, klev
739       DO i = 1, klon
740        omega(i,k) = RG*flxmw(i,k) / cell_area(i)
741       END DO
742      END DO
743
744c Ajouter le geopotentiel du sol:
745c
746      DO k = 1, klev
747      DO i = 1, klon
748         zphi(i,k) = pphi(i,k) + pphis(i)
749      ENDDO
750      ENDDO
751
752c     call WriteField_phy('physiq_pphi',pphi,klev)
753c     call WriteField_phy('physiq_pphis',pphis,1)
754
755c   calcul du geopotentiel aux niveaux intercouches
756c   ponderation des altitudes au niveau des couches en dp/p
757
758      DO l=1,klev
759         DO i=1,klon
760c           zzlay(i,l)=zphi(i,l)/RG
761c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
762            zzlay(i,l)=RG*RA*RA/(RG*RA-zphi(i,l))-RA
763         ENDDO
764      ENDDO
765      DO i=1,klon
766c        zzlev(i,1)=0.
767c CORRECTION 13/01/2011 
768c (correspond a la position de la surface en ce point vs RA)
769         zzlev(i,1)=pphis(i)/RG
770      ENDDO
771      DO l=2,klev
772         DO i=1,klon
773            z1=(pplay(i,l-1)+paprs(i,l))/(pplay(i,l-1)-paprs(i,l))
774            z2=(paprs(i,l)  +pplay(i,l))/(paprs(i,l)  -pplay(i,l))
775            zzlev(i,l)=(z1*zzlay(i,l-1)+z2*zzlay(i,l))/(z1+z2)
776         ENDDO
777      ENDDO
778      DO i=1,klon
779         zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
780      ENDDO
781
782! zonal averages needed
783      if (moyzon_ch.or.moyzon_mu) then
784
785c      zzlaybar(1,:)=(zphibar(1,:)+zphisbar(1))/RG
786c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
787       zzlaybar(1,:)=RG*RA*RA/(RG*RA-(zphibar(1,:)+zphisbar(1)))-RA
788       zzlevbar(1,1)=zphisbar(1)/RG
789       DO l=2,klev
790            z1=(zplaybar(1,l-1)+zplevbar(1,l))/
791     .            (zplevbar(1,l-1)-zplevbar(1,l))
792            z2=(zplevbar(1,l)  +zplaybar(1,l))/
793     .            (zplevbar(1,l)  -zplaybar(1,l))
794            zzlevbar(1,l)=(z1*zzlaybar(1,l-1)+z2*zzlaybar(1,l))/(z1+z2)
795       ENDDO
796       zzlevbar(1,klev+1)=zzlaybar(1,klev)+
797     .            (zzlaybar(1,klev)-zzlevbar(1,klev))
798
799       DO i=2,klon
800        if (latitude(i).ne.latitude(i-1)) then
801         DO l=1,klev
802c         zzlaybar(i,l)=(zphibar(i,l)+zphisbar(i))/RG
803c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
804          zzlaybar(i,l)=RG*RA*RA/(RG*RA-(zphibar(i,l)+zphisbar(i)))-RA
805         ENDDO
806         zzlevbar(i,1)=zphisbar(i)/RG
807         DO l=2,klev
808            z1=(zplaybar(i,l-1)+zplevbar(i,l))/
809     .            (zplevbar(i,l-1)-zplevbar(i,l))
810            z2=(zplevbar(i,l)  +zplaybar(i,l))/
811     .            (zplevbar(i,l)  -zplaybar(i,l))
812            zzlevbar(i,l)=(z1*zzlaybar(i,l-1)+z2*zzlaybar(i,l))/(z1+z2)
813         ENDDO
814         zzlevbar(i,klev+1)=zzlaybar(i,klev)+
815     .              (zzlaybar(i,klev)-zzlevbar(i,klev))
816        else
817         zzlaybar(i,:)=zzlaybar(i-1,:)
818         zzlevbar(i,:)=zzlevbar(i-1,:)
819        endif
820       ENDDO
821
822      endif  ! moyzon
823
824c     call WriteField_phy('physiq_zphi',zphi,klev)
825c     call WriteField_phy('physiq_zzlay',zzlay,klev)
826c     call WriteField_phy('physiq_zzlev',zzlev,klev+1)
827c- - - - - - - - - - - - - - - -
828c DIAGNOSTIQUE GRILLE VERTICALE
829c- - - - - - - - - - - - - - - -
830c     print*,"DIAGNOSTIQUE GRILLE VERTICALE"
831c     i=klon/2
832c     print*,"Niveau  Pression  Altitude    (lev puis lay)"
833c     do l=1,klev
834c      print*,l,paprs(i,l),zzlev(i,l)
835c      print*,l,pplay(i,l),zzlay(i,l)
836c     enddo
837c     print*,klev+1,paprs(i,klev+1),zzlev(i,klev+1)
838c     stop
839
840c====================================================================
841c
842c Verifier les temperatures
843c
844      CALL hgardfou(t_seri,ftsol,'debutphy')
845c====================================================================
846c
847c Incrementer le compteur de la physique
848c
849      itap   = itap + 1
850
851c====================================================================
852c
853c Epaisseurs couches
854
855      DO k = 1, klev
856      DO i = 1, klon
857         delp(i,k) = paprs(i,k)-paprs(i,k+1)
858      ENDDO
859      ENDDO
860
861c====================================================================
862c Orbite et eclairement
863c====================================================================
864
865c Pour TITAN:
866c  calcul de la longitude solaire
867          CALL solarlong(rjourvrai+gmtime,zls)
868          zlsdeg = zls*180./RPI      ! zls est en radians !!
869          print*,'Ls',zlsdeg
870
871      CALL orbite(zls,dist,pdecli)
872      IF (debut) zlsm1=zls
873
874c dans zenang, Ls en degres ; dans mucorr, Ls en radians
875      call mucorr(klon,zls,latitude_deg,rmu0bar,fractbar)
876      IF (cycle_diurne) THEN
877        zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
878        CALL zenang(zlsdeg,gmtime,zdtime,latitude_deg,longitude_deg,
879     &              rmu0,fract)
880      ELSE
881        rmu0  = rmu0bar
882        fract = fractbar
883      ENDIF
884     
885c====================================================================
886c Appeler la diffusion verticale (programme de couche limite)
887c====================================================================
888
889c-------------------------------
890c TEST: on ne tient pas compte des calculs de clmain mais on force
891c l'equilibre radiatif du sol
892      if (1.eq.0) then
893              if (debut) then
894                print*,"ATTENTION, CLMAIN SHUNTEE..."
895              endif
896
897      DO i = 1, klon
898         sens(i) = 0.0e0 ! flux de chaleur sensible au sol
899         fder(i) = 0.0e0
900         dlw(i)  = 0.0e0
901      ENDDO
902
903c Incrementer la temperature du sol
904c
905      DO i = 1, klon
906         d_ts(i)  = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
907         ftsol(i) = ftsol(i) + d_ts(i)
908         do j=1,nsoilmx
909           ftsoil(i,j)=ftsol(i)
910         enddo
911      ENDDO
912
913c-------------------------------
914      else
915c-------------------------------
916
917      fder = dlw
918
919c     print*,"radsol avant clmain=",radsol(klon/2)
920c     print*,"solsw avant clmain=",solsw(klon/2)
921c     print*,"sollw avant clmain=",sollw(klon/2)
922
923! ADAPTATION GCM POUR CP(T)
924
925      CALL clmain(dtime,itap,
926     e            t_seri,u_seri,v_seri,
927     e            rmu0,
928     e            ftsol,
929     $            ftsoil,
930     $            paprs,pplay,ppk,radsol,falbe,
931     e            solsw, sollw, sollwdown, fder,
932     e            longitude_deg, latitude_deg, dx, dy,   
933     e            debut, lafin,
934     s            d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
935     s            fluxt,fluxu,fluxv,cdragh,cdragm,
936     s            dsens,
937     s            ycoefh,yu1,yv1)
938
939c     print*,"radsol apres clmain=",radsol(klon/2)
940c     print*,"solsw apres clmain=",solsw(klon/2)
941c     print*,"sollw apres clmain=",sollw(klon/2)
942
943CXXX Incrementation des flux
944      DO i = 1, klon
945         sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
946         fder(i) = dlw(i) + dsens(i)
947      ENDDO
948CXXX
949
950      DO k = 1, klev
951      DO i = 1, klon
952         t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
953         d_t_vdf(i,k)= d_t_vdf(i,k)/dtime          ! K/s
954         u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
955         d_u_vdf(i,k)= d_u_vdf(i,k)/dtime          ! (m/s)/s
956         v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
957         d_v_vdf(i,k)= d_v_vdf(i,k)/dtime          ! (m/s)/s
958      ENDDO
959      ENDDO
960
961c     call WriteField_phy('physiq_dtvdf',d_t_vdf,klev)
962c     call WriteField_phy('physiq_duvdf',d_u_vdf,klev)
963c     call WriteField_phy('physiq_dvvdf',d_v_vdf,klev)
964
965C TRACEURS
966
967      d_tr_vdf = 0.
968      if (iflag_trac.eq.1) then
969         DO iq=1, nqmax
970             CALL cltrac(dtime,ycoefh,t_seri,
971     s               tr_seri(1,1,iq),source,
972     e               paprs, pplay,delp,
973     s               d_tr_vdf(1,1,iq))
974             tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
975             d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime          ! /s
976         ENDDO
977      endif
978
979      IF (if_ebil.ge.2) THEN
980        ztit='after clmain'
981        CALL diagetpq(cell_area,ztit,ip_ebil,2,1,dtime
982     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
983     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
984         call diagphy(cell_area,ztit,ip_ebil
985     e      , zero_v, zero_v, zero_v, zero_v, sens
986     e      , zero_v, zero_v, zero_v, ztsol
987     e      , d_h_vcol, d_qt, d_ec
988     s      , fs_bound, fq_bound )
989      END IF
990C
991c
992c Incrementer la temperature du sol
993c
994c     print*,'Tsol avant clmain:',ftsol(1)
995      DO i = 1, klon
996         ftsol(i) = ftsol(i) + d_ts(i)
997      ENDDO
998c     print*,'DTsol apres clmain:',d_ts(klon/2)
999c     print*,'Tsol apres clmain:',ftsol(1)
1000
1001c Calculer la derive du flux infrarouge
1002c
1003      DO i = 1, klon
1004            dlw(i) = - 4.0*emis*RSIGMA*ftsol(i)**3
1005      ENDDO
1006
1007c-------------------------------
1008      endif  ! fin du TEST
1009
1010c
1011c Appeler l'ajustement sec
1012c
1013c===================================================================
1014c Convection seche
1015c===================================================================
1016c
1017      d_t_ajs(:,:)=0.
1018      d_u_ajs(:,:)=0.
1019      d_v_ajs(:,:)=0.
1020      d_tr_ajs(:,:,:)=0.
1021c
1022      IF(prt_level>9)WRITE(lunout,*)
1023     .    'AVANT LA CONVECTION SECHE , iflag_ajs='
1024     s   ,iflag_ajs
1025
1026      if(iflag_ajs.eq.0) then
1027c  Rien
1028c  ====
1029         IF(prt_level>9)WRITE(lunout,*)'pas de convection'
1030
1031      else if(iflag_ajs.eq.1) then
1032
1033c  Ajustement sec
1034c  ==============
1035         IF(prt_level>9)WRITE(lunout,*)'ajsec'
1036
1037! ADAPTATION GCM POUR CP(T)
1038         CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
1039     .              tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
1040
1041! ADAPTATION GCM POUR CP(T)
1042         do i=1,klon
1043          flux_ajs(i,1) = 0.0
1044          do j=2,klev
1045            flux_ajs(i,j) = flux_ajs(i,j-1)
1046     .        + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime
1047     .                                 *delp(i,j-1)
1048          enddo
1049         enddo
1050         
1051         t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
1052         d_t_ajs(:,:)= d_t_ajs(:,:)/dtime          ! K/s
1053         u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
1054         d_u_ajs(:,:)= d_u_ajs(:,:)/dtime          ! (m/s)/s
1055         v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
1056         d_v_ajs(:,:)= d_v_ajs(:,:)/dtime          ! (m/s)/s
1057      if (iflag_trac.eq.1) then
1058           tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
1059           d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime  ! /s
1060      endif
1061
1062c     call WriteField_phy('physiq_dtajs',d_t_ajs,klev)
1063c     call WriteField_phy('physiq_duajs',d_u_ajs,klev)
1064c     call WriteField_phy('physiq_dvajs',d_v_ajs,klev)
1065
1066      endif
1067c
1068      IF (if_ebil.ge.2) THEN
1069        ztit='after dry_adjust'
1070        CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
1071     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1072     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1073        call diagphy(cell_area,ztit,ip_ebil
1074     e      , zero_v, zero_v, zero_v, zero_v, sens
1075     e      , zero_v, zero_v, zero_v, ztsol
1076     e      , d_h_vcol, d_qt, d_ec
1077     s      , fs_bound, fq_bound )
1078      END IF
1079
1080c====================================================================
1081c   MICROPHYSIQUE ET CHIMIE
1082c====================================================================
1083
1084      d_tr_mph(:,:,:)=0.
1085      d_tr_kim(:,:,:)=0.
1086     
1087c on recupere tr_seri inchange, d_tr_micro, d_tr_chim, tous les trois sur nqmax
1088c on recupere aussi qaer pour le mettre dans les sorties
1089c  si microfi=1, sortie de qaer(1:nmicro)
1090c  si nmicro != nqmax et si chimi, sortie de tr_seri(nmicro+1:nqmax)
1091
1092c faire un test comme pour rayonnement, avec chimi en plus comme flag,
1093c pour voir si chimie appelee -> bouleen, qui passe dans phytrac.
1094c faut aussi le pas de temps chimique: dtimechim, a passer..
1095
1096      appel_chim = 0
1097      IF (MOD(itapchim,chimpas).EQ.0) THEN
1098c             print*,'CHIMIE ',
1099c    $             ' (itapchim=',itapchim,'/chimpas=',chimpas,')'
1100       appel_chim = 1
1101       itapchim = 0
1102      ENDIF
1103      itapchim = itapchim + 1
1104
1105      if (iflag_trac.eq.1) then
1106c        call WriteField_phy('physiq_qaer01',
1107c    .                          qaer(:,:,1),klev)
1108c        call WriteField_phy('physiq_qaer10',
1109c    .                          qaer(:,:,10),klev)
1110c        call WriteField_phy('physiq_tr_seri01',
1111c    .                          tr_seri(:,:,1),klev)
1112c        call WriteField_phy('physiq_tr_seri10',
1113c    .                          tr_seri(:,:,10),klev)
1114
1115c         call begintime(tt0)
1116c in phytrac call, mu0 and fract are only used in brume
1117c so we need to pass either rmu0 ou rmu0bar depending on
1118c moyzon_mu
1119       if (moyzon_mu) then
1120         call phytrac (debut,lafin,
1121     .                 nqmax,nmicro,dtime,appel_chim,dtimechim,
1122     .                 paprs,pplay,delp,t,rmu0bar,fractbar,pdecli,zls,
1123     .                 yu1,yv1,zzlev,zzlay,ftsol,
1124     .                 tr_seri,qaer,d_tr_mph,d_tr_kim,
1125     .                 fclat,reservoir)
1126       else
1127         call phytrac (debut,lafin,
1128     .                 nqmax,nmicro,dtime,appel_chim,dtimechim,
1129     .                 paprs,pplay,delp,t,rmu0,fract,pdecli,zls,
1130     .                 yu1,yv1,zzlev,zzlay,ftsol,
1131     .                 tr_seri,qaer,d_tr_mph,d_tr_kim,
1132     .                 fclat,reservoir)
1133       endif
1134
1135c         call endtime(tt0,tt1)
1136c         ttphytra=ttphytra+tt1
1137
1138c ----- ICI on ajuste radsol en tenant compte du flux de chaleur latente
1139c       d'evaporation du reservoir.
1140c       NOTE : c'est pas tres elegant mais ca permet d'eviter d'aller
1141c              toucher a clmain.
1142        if (clouds.eq.1) then
1143          radsol(:) = radsol(:)+fclat(:)    !test pas de flx de chaleur latente
1144        endif
1145
1146        if (microfi.ge.1) then
1147         tr_seri(:,:,1:nmicro) = tr_seri(:,:,1:nmicro)
1148     .                        + d_tr_mph(:,:,1:nmicro)*dtime
1149c        call WriteField_phy('physiq_d_tr_mph01',
1150c    .                          d_tr_mph(:,:,1),klev)
1151c        call WriteField_phy('physiq_d_tr_mph10',
1152c    .                          d_tr_mph(:,:,10),klev)
1153        endif
1154c       PAS ELEGANT mais je n'ai pas trouve d'autres solutions :
1155c       Il semblerait qu'il y ait un probleme lorsque les tendances de traceurs
1156c       retourne des traceurs nuls et il y a parfois des valeurs negatives qui trainent.
1157c       Pour ne diffuser le probleme, on force les valeurs negatives a ZERO.
1158        DO iq=1,nmicro
1159          DO i=1,klon
1160            DO l=1,klev
1161              if (tr_seri(i,l,iq).lt.0.) then
1162                 tr_seri(i,l,iq) = 0.
1163              endif
1164            ENDDO
1165          ENDDO
1166        ENDDO
1167
1168c condensation:
1169c       NE PAS OUBLIER LA CONDENSATION DES NUAGES !!!!
1170        if ((clouds.eq.1.or.(chimi)).and.nqmax.gt.nmicro) then
1171          tr_seri(:,:,nmicro+1:nqmax) = tr_seri(:,:,nmicro+1:nqmax)
1172     .                         + d_tr_mph(:,:,nmicro+1:nqmax)*dtime
1173        endif
1174
1175c chimie:
1176        if ((chimi).and.(nqmax.gt.nmicro)) then
1177         tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_kim(:,:,:)*dtime
1178        endif
1179
1180      endif  !iflag_trac=1
1181
1182c       ch4=0.
1183c       do l=1,llm
1184c         ch4(1,l) = tr_seri(1,l,ich4)
1185c         do j=2,jjm
1186c           ig0=1+(j-2)*iim
1187c           do i=1,iim
1188c             ch4(j,l)= ch4(j,l)  + tr_seri(ig0+i,l,ich4)/iim
1189c           enddo
1190c         enddo
1191c         ch4(jjm+1,l) = tr_seri(klon,l,ich4)
1192c       enddo
1193c       do j=1,jjm+1
1194c         write(501,*) j,ch4(j,1)
1195c       enddo
1196c       do l=1,llm
1197c         write(502,'(I3,49(ES24.17,1X))') l,
1198c     &   (ch4(j,l),j=1,jjm+1)
1199c       enddo
1200c       write(501,*) ""
1201c       write(502,*) ""
1202
1203c------------------
1204c test condensation
1205c     do i=1,nqmax
1206c       if(tname(i).eq."HCN") then
1207c          print*,"HCN="
1208c          do k=1,klev
1209c           print*,k,tr_seri(klon/2,k,i),d_tr_mph(klon/2,k,i)*dtime
1210c    v      ,d_tr_kim(klon/2,k,i)*dtime
1211c          enddo
1212c          stop
1213c       endif
1214c     enddo
1215c------------------
1216
1217c====================================================================
1218c RAYONNEMENT
1219c====================================================================
1220
1221      IF (MOD(itaprad,radpas).EQ.0) THEN
1222c             print*,'RAYONNEMENT ',
1223c    $             ' (itaprad=',itaprad,'/radpas=',radpas,')'
1224
1225c ATTENTION, (klon/2) ne marche pas toujours............
1226c     print*,"radsol avant radlwsw=",radsol(klon/2)
1227c     print*,"solsw avant radlwsw=",solsw(klon/2)
1228c     print*,"sollw avant radlwsw=",sollw(klon/2)
1229c     print*,"avant radlwsw"
1230
1231c   ----------------
1232c   Calcul du rayon moyen des gouttes et des fractions volumique pour le TR
1233c  ----------------
1234      IF (clouds.eq.1) THEN
1235        DO i=1,klon
1236          DO j=1,klev
1237            rmcbar(i,j)=rmcbar(i,j)/MAX(REAL(ncount(i,j)),1.)
1238            xfbar(i,j,:)=xfbar(i,j,:)/MAX(REAL(ncount(i,j)),1.)
1239          ENDDO
1240        ENDDO
1241      ENDIF
1242     
1243c      call begintime(tt0)
1244      CALL radlwsw
1245     e            (dist, rmu0, fract, zzlev,
1246     e             paprs, pplay,ftsol, t_seri, nqmax, nmicro,
1247     c             tr_seri, qaer)
1248c     print*,"apres radlwsw"
1249
1250c      call endtime(tt0,tt1)
1251c      ttrad=ttrad+tt1
1252
1253c     print*,"apres radlwsw"
1254c     mise a zero du rayon moyen des gouttes et des fractions volumique
1255      IF (clouds.eq.1) THEN
1256        rmcbar(:,:)  = 0.
1257        xfbar(:,:,:) = 0.
1258        ncount(:,:)  = 0
1259      ENDIF
1260
1261      itaprad = 0
1262      DO k = 1, klev
1263       DO i = 1, klon
1264         dtrad(i,k) = heat(i,k)-cool(i,k)     !K/s
1265       ENDDO
1266      ENDDO
1267
1268c     call WriteField_phy('physiq_heat',heat,klev)
1269c     call WriteField_phy('physiq_cool',cool,klev)
1270
1271      ENDIF
1272      itaprad = itaprad + 1
1273c====================================================================
1274c
1275c Ajouter la tendance des rayonnements (tous les pas)
1276c
1277      DO k = 1, klev
1278      DO i = 1, klon
1279         t_seri(i,k) = t_seri(i,k) + dtrad(i,k) * dtime
1280      ENDDO
1281      ENDDO
1282 
1283      IF (if_ebil.ge.2) THEN
1284        ztit='after rad'
1285        CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
1286     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1287     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1288        call diagphy(cell_area,ztit,ip_ebil
1289     e      , topsw, toplw, solsw, sollw, zero_v
1290     e      , zero_v, zero_v, zero_v, ztsol
1291     e      , d_h_vcol, d_qt, d_ec
1292     s      , fs_bound, fq_bound )
1293      END IF
1294c
1295
1296c====================================================================
1297c   Calcul  des gravity waves  FLOTT
1298c====================================================================
1299c
1300      if (ok_orodr.or.ok_gw_nonoro) then
1301c  CALCUL DE N2
1302       do i=1,klon
1303        do k=2,klev
1304          ztlev(i,k)  = (t_seri(i,k)+t_seri(i,k-1))/2.
1305          zpklev(i,k) = sqrt(ppk(i,k)*ppk(i,k-1))
1306        enddo
1307       enddo
1308       call t2tpot(klon*klev,ztlev, ztetalev,zpklev)
1309       call t2tpot(klon*klev,t_seri,ztetalay,ppk)
1310       do i=1,klon
1311        do k=2,klev
1312          zdtetalev(i,k) = ztetalay(i,k)-ztetalay(i,k-1)
1313          zdzlev(i,k)    = (zphi(i,k)-zphi(i,k-1))/RG
1314          zn2(i,k) = RG*zdtetalev(i,k)/(ztetalev(i,k)*zdzlev(i,k))
1315          zn2(i,k) = max(zn2(i,k),1.e-12)  ! securite
1316        enddo
1317        zn2(i,1) = 1.e-12  ! securite
1318       enddo
1319
1320      endif
1321     
1322c ----------------------------ORODRAG
1323      IF (ok_orodr) THEN
1324c
1325c  selection des points pour lesquels le shema est actif:
1326        igwd=0
1327        DO i=1,klon
1328        itest(i)=0
1329c        IF ((zstd(i).gt.10.0)) THEN
1330        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1331          itest(i)=1
1332          igwd=igwd+1
1333          idx(igwd)=i
1334        ENDIF
1335        ENDDO
1336c        igwdim=MAX(1,igwd)
1337c
1338c A ADAPTER POUR VENUS!!!
1339        CALL drag_noro(klon,klev,dtime,paprs,pplay,pphi,zn2,
1340     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
1341     e                   igwd,idx,itest,
1342     e                   t_seri, u_seri, v_seri,
1343     s                   zulow, zvlow, zustrdr, zvstrdr,
1344     s                   d_t_oro, d_u_oro, d_v_oro)
1345
1346c       print*,"d_u_oro=",d_u_oro(klon/2,:)
1347c  ajout des tendances
1348           t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
1349           d_t_oro(:,:)= d_t_oro(:,:)/dtime          ! K/s
1350           u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
1351           d_u_oro(:,:)= d_u_oro(:,:)/dtime          ! (m/s)/s
1352           v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
1353           d_v_oro(:,:)= d_v_oro(:,:)/dtime          ! (m/s)/s
1354c   
1355      ELSE
1356         d_t_oro = 0.
1357         d_u_oro = 0.
1358         d_v_oro = 0.
1359         zustrdr = 0.
1360         zvstrdr = 0.
1361c
1362      ENDIF ! fin de test sur ok_orodr
1363c
1364c ----------------------------OROLIFT
1365      IF (ok_orolf) THEN
1366       print*,"ok_orolf NOT IMPLEMENTED !"
1367       stop
1368c
1369c  selection des points pour lesquels le shema est actif:
1370        igwd=0
1371        DO i=1,klon
1372        itest(i)=0
1373        IF ((zpic(i)-zmea(i)).GT.100.) THEN
1374          itest(i)=1
1375          igwd=igwd+1
1376          idx(igwd)=i
1377        ENDIF
1378        ENDDO
1379c        igwdim=MAX(1,igwd)
1380c
1381c A ADAPTER POUR VENUS ET TITAN!!!
1382c            CALL lift_noro(klon,klev,dtime,paprs,pplay,
1383c     e                   latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval,
1384c     e                   igwd,idx,itest,
1385c     e                   t_seri, u_seri, v_seri,
1386c     s                   zulow, zvlow, zustrli, zvstrli,
1387c     s                   d_t_lif, d_u_lif, d_v_lif               )
1388
1389c
1390c  ajout des tendances
1391           t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
1392           d_t_lif(:,:)= d_t_lif(:,:)/dtime          ! K/s
1393           u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
1394           d_u_lif(:,:)= d_u_lif(:,:)/dtime          ! (m/s)/s
1395           v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
1396           d_v_lif(:,:)= d_v_lif(:,:)/dtime          ! (m/s)/s
1397c
1398      ELSE
1399         d_t_lif = 0.
1400         d_u_lif = 0.
1401         d_v_lif = 0.
1402         zustrli = 0.
1403         zvstrli = 0.
1404c
1405      ENDIF ! fin de test sur ok_orolf
1406
1407c ---------------------------- NON-ORO GRAVITY WAVES
1408       IF(ok_gw_nonoro) then
1409
1410        abort_message="Option non developpee pour Titan"
1411        call abort_gcm(modname,abort_message,1)
1412c A FAIRE POUR TITAN
1413c      call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
1414c     e               t_seri, u_seri, v_seri,
1415c     o               zustrhi,zvstrhi,
1416c     o               d_t_hin, d_u_hin, d_v_hin)
1417
1418c  ajout des tendances
1419
1420c         t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
1421c         d_t_hin(:,:)= d_t_hin(:,:)/dtime          ! K/s
1422c         u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
1423c         d_u_hin(:,:)= d_u_hin(:,:)/dtime          ! (m/s)/s
1424c         v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
1425c         d_v_hin(:,:)= d_v_hin(:,:)/dtime          ! (m/s)/s
1426
1427      ELSE
1428         d_t_hin = 0.
1429         d_u_hin = 0.
1430         d_v_hin = 0.
1431         zustrhi = 0.
1432         zvstrhi = 0.
1433
1434      ENDIF ! fin de test sur ok_gw_nonoro
1435
1436c====================================================================
1437c Transport de ballons
1438c====================================================================
1439      if (ballons.eq.1) then
1440         CALL ballon(30,pdtphys,rjourvrai,gmtime,
1441     &               latitude_deg,longitude_deg,
1442c    C               t,pplay,u,v,pphi)   ! alt above surface (smoothed for GCM)
1443     C               t,pplay,u,v,zphi)   ! alt above planet average radius
1444      endif !ballons
1445
1446c====================================================================
1447c Bilan de mmt angulaire
1448c====================================================================
1449      if (bilansmc.eq.1) then
1450CMODDEB FLOTT
1451C  CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
1452C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
1453
1454      DO i = 1, klon
1455        zustrph(i)=0.
1456        zvstrph(i)=0.
1457        zustrcl(i)=0.
1458        zvstrcl(i)=0.
1459      ENDDO
1460      DO k = 1, klev
1461      DO i = 1, klon
1462       zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
1463     c            (paprs(i,k)-paprs(i,k+1))/rg
1464       zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
1465     c            (paprs(i,k)-paprs(i,k+1))/rg
1466       zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
1467     c            (paprs(i,k)-paprs(i,k+1))/rg
1468       zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
1469     c            (paprs(i,k)-paprs(i,k+1))/rg
1470      ENDDO
1471      ENDDO
1472
1473      CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
1474     C               ra,rg,romega,
1475     C               latitude_deg,longitude_deg,pphis,
1476     C               zustrdr,zustrli,zustrcl,
1477     C               zvstrdr,zvstrli,zvstrcl,
1478     C               paprs,u,v)
1479                     
1480CCMODFIN FLOTT
1481      endif !bilansmc
1482
1483c====================================================================
1484c====================================================================
1485c Calculer le transport de l'eau et de l'energie (diagnostique)
1486c
1487c  A REVOIR POUR VENUS ET TITAN...
1488c
1489c     CALL transp (paprs,ftsol,
1490c    e                   t_seri, q_seri, u_seri, v_seri, zphi,
1491c    s                   ve, vq, ue, uq)
1492c
1493c====================================================================
1494c+jld ec_conser
1495      DO k = 1, klev
1496      DO i = 1, klon
1497        d_t_ec(i,k)=0.5/cpdet(t_seri(i,k))
1498     $      *(u(i,k)**2+v(i,k)**2-u_seri(i,k)**2-v_seri(i,k)**2)
1499        t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
1500        d_t_ec(i,k) = d_t_ec(i,k)/dtime
1501       END DO
1502      END DO
1503         do i=1,klon
1504          flux_ec(i,1) = 0.0
1505          do j=2,klev
1506            flux_ec(i,j) = flux_ec(i,j-1)
1507     . +cpdet(t_seri(i,j-1))/RG*d_t_ec(i,j-1)*delp(i,j-1)
1508          enddo
1509         enddo
1510         
1511c-jld ec_conser
1512c====================================================================
1513      IF (if_ebil.ge.1) THEN
1514        ztit='after physic'
1515        CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
1516     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1517     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1518C     Comme les tendances de la physique sont ajoute dans la dynamique,
1519C     on devrait avoir que la variation d'entalpie par la dynamique
1520C     est egale a la variation de la physique au pas de temps precedent.
1521C     Donc la somme de ces 2 variations devrait etre nulle.
1522        call diagphy(cell_area,ztit,ip_ebil
1523     e      , topsw, toplw, solsw, sollw, sens
1524     e      , zero_v, zero_v, zero_v, ztsol
1525     e      , d_h_vcol, d_qt, d_ec
1526     s      , fs_bound, fq_bound )
1527C
1528      d_h_vcol_phy=d_h_vcol
1529C
1530      END IF
1531C
1532c=======================================================================
1533c   SORTIES
1534c=======================================================================
1535
1536c Convertir les incrementations en tendances
1537c
1538      DO k = 1, klev
1539      DO i = 1, klon
1540         d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
1541         d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
1542         d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
1543      ENDDO
1544      ENDDO
1545c
1546      DO iq = 1, nqmax
1547      DO  k = 1, klev
1548      DO  i = 1, klon
1549         d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
1550      ENDDO
1551      ENDDO
1552      ENDDO
1553     
1554c------------------------
1555c Calcul moment cinetique
1556c------------------------
1557c TEST...
1558c     mangtot = 0.0
1559c     DO k = 1, klev
1560c     DO i = 1, klon
1561c       mang(i,k) = RA*cos(latitude(i))
1562c    .     *(u_seri(i,k)+RA*cos(latitude(i))*ROMEGA)
1563c    .     *cell_area(i)*(paprs(i,k)-paprs(i,k+1))/RG
1564c       mangtot=mangtot+mang(i,k)
1565c     ENDDO
1566c     ENDDO
1567c     print*,"Moment cinetique total = ",mangtot
1568c
1569c------------------------
1570c
1571c Sauvegarder les valeurs de t et u a la fin de la physique:
1572c
1573      DO k = 1, klev
1574      DO i = 1, klon
1575         u_ancien(i,k) = u_seri(i,k)
1576         t_ancien(i,k) = t_seri(i,k)
1577      ENDDO
1578      ENDDO
1579c
1580c=============================================================
1581c   Ecriture des sorties
1582c=============================================================
1583     
1584#ifdef CPP_IOIPSL
1585
1586#ifdef histday
1587#include "write_histday.h"
1588#endif
1589
1590#ifdef histmth
1591#include "write_histmth.h"
1592#endif
1593
1594#ifdef histins
1595#include "write_histins.h"
1596#endif
1597
1598#endif
1599
1600c====================================================================
1601c Si c'est la fin, il faut conserver l'etat de redemarrage
1602c====================================================================
1603c
1604      IF (lafin) THEN
1605         itau_phy = itau_phy + itap
1606         lsinit   = zlsdeg
1607         CALL phyredem ("restartphy.nc")
1608     
1609c--------------FLOTT
1610CMODEB LOTT
1611C  FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
1612C  DU BILAN DE MOMENT ANGULAIRE.
1613      if (bilansmc.eq.1) then
1614        write(*,*)'Fermeture de aaam_bud.out (FL Vous parle)'
1615        close(27)                                     
1616        close(28)                                     
1617      endif !bilansmc
1618CMODFIN
1619c-------------
1620c--------------SLEBONNOIS
1621C  FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
1622C  DES BALLONS
1623      if (ballons.eq.1) then
1624        write(*,*)'Fermeture des ballons*.out'
1625        close(30)                                     
1626        close(31)                                     
1627        close(32)                                     
1628        close(33)                                     
1629        close(34)                                     
1630      endif !ballons
1631c-------------
1632      ENDIF
1633     
1634      END SUBROUTINE physiq
1635
1636      END MODULE physiq_mod
1637
Note: See TracBrowser for help on using the repository browser.