source: LMDZ6/branches/LMDZ-tracers/libf/phymar/physiq_mod.F90 @ 3999

Last change on this file since 3999 was 2599, checked in by Ehouarn Millour, 8 years ago

Update for phymar which should have been done with revision 2315.
Introduce module vertical_layers_mod to store information on the
vertical grid to be used within the physics (rather than including
the dynamics-related comvert.h).
EM

File size: 67.2 KB
Line 
1! $Id: physiq.F 1565 2011-08-31 12:53:29Z jghattas $
2!#define IO_DEBUG
3MODULE physiq_mod
4
5IMPLICIT NONE
6
7CONTAINS
8
9!========================================================================================================================
10! Adaptation de physiq.F90 de phydev pour mettre en oeuvre une interface dynamique de LMDZ/physique de MAR.
11! Martin Ménégoz (Février 2014) martin.menegoz@lgge.obs.ujf-grenoble.fr
12! Contact LGGE: gallee@lgge.obs.ujf-grenoble.fr (modèle MAR) ; contact LMD: Frederic.hourdin@lmd.jussieu.fr (modèle LMDZ)
13!========================================================================================================================
14
15      SUBROUTINE physiq (nlon,nlev, &
16     &            debut,lafin,jD_cur, jH_cur,pdtphys, &
17     &            paprs,pplay,pphi,pphis,ppresnivs, &
18     &            u,v,rot,t,qx, &
19     &            flxmass_w, &
20     &            d_u, d_v, d_t, d_qx, d_ps &
21     &            , dudyn )
22
23      USE dimphy, only : klon,klev,klevp1
24      USE infotrac_phy, only : nqtot
25      USE geometry_mod, only : latitude,longitude,cell_area
26      !USE comcstphy, only : rg
27      USE iophy, only : histbeg_phy,histwrite_phy
28      USE ioipsl, only : getin,histvert,histdef,histend,ymds2ju,ju2ymds,histsync
29      USE mod_phys_lmdz_para, only : jj_nb
30      USE phys_state_var_mod, only : phys_state_var_init
31      USE vertical_layers_mod, ONLY: ap, bp
32! Modules LMDZ utilisés pour LMDZ-MAR:
33      USE control_mod     
34
35! Modules utilisé pour MAR dans DYN_to_PHY_MAR.f90
36      use Mod_Real
37      use Mod_PHY____dat
38      use Mod_PHY____grd
39      use Mod_PHY_S0_grd
40      use Mod_SISVAT_grd
41      use Mod_SISVAT_dat
42      use Mod_SISVAT_gpt ! surface heat flxes in outputs
43      use Mod_PHY_DY_kkl
44      use Mod_PHY_CM_kkl ! pour écrire les snowfall et rainfall dans les outputs
45      use Mod_PHY_CP_kkl ! pour écrire les snowfall et rainfall dans les outputs
46
47#ifdef CPP_XIOS
48      USE wxios
49#endif
50
51      IMPLICIT none
52
53#include "dimensions.h"
54! Include LMDZ utilisés ici pour MAR-LMDZ:
55#include "YOMCST.h"
56
57!======================================================================
58! Arguments du moniteur general de la physique dans le modele LMDZ
59!======================================================================
60!
61!  Arguments:
62!
63! nlon----input-I-nombre de points horizontaux
64! nlev----input-I-nombre de couches verticales, doit etre egale a klev
65! debut---input-L-variable logique indiquant le premier passage
66! lafin---input-L-variable logique indiquant le dernier passage
67! jD_cur       -R-jour courant a l'appel de la physique (jour julien)
68! jH_cur       -R-heure courante a l'appel de la physique (jour julien)
69! pdtphys-input-R-pas d'integration pour la physique (seconde)
70! paprs---input-R-pression pour chaque inter-couche (en Pa)
71! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
72! pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
73! pphis---input-R-geopotentiel du sol
74! ppresnivs-input_R_pressions approximat. des milieux couches ( en PA)
75! u-------input-R-vitesse dans la direction X (de O a E) en m/s
76! v-------input-R-vitesse Y (de S a N) en m/s
77! t-------input-R-temperature (K)
78! qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
79! d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
80! d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
81! flxmass_w -input-R- flux de masse verticale
82! d_u-----output-R-tendance physique de "u" (m/s/s)
83! d_v-----output-R-tendance physique de "v" (m/s/s)
84! d_t-----output-R-tendance physique de "t" (K/s)
85! d_qx----output-R-tendance physique de "qx" (kg/kg/s)
86! d_ps----output-R-tendance physique de la pression au sol
87!IM
88
89!======================================================================
90! Variables ajoutés pour l'utilisation de la physique de MAR
91!======================================================================
92
93      integer,parameter :: jjmp1=jjm+1-1/jjm
94      integer,parameter :: iip1=iim+1
95
96! Indices des traceurs (qx)
97! NB: Tous les traceurs de MAR (vapeur d'eau et hydrométéores, solides et liquides) sont transportés via la variables qx (tendances: d_qx)
98! Il doivent être déclarés dans un fichier traceur.def avant de lancer la simulation
99
100      INTEGER, PARAMETER :: ivap  = 1  ! indice de traceurs pour vapeur d'eau
101      INTEGER, PARAMETER :: iliq  = 2  ! indice de traceurs pour eau liquide
102      INTEGER, PARAMETER :: iice  = 3  ! indice de traceurs pour eau solide
103      INTEGER, PARAMETER :: icin  = 4  ! indice de traceurs pour CIN
104      INTEGER, PARAMETER :: isnow = 5  ! indice de traceurs pour flocons de neige
105      INTEGER, PARAMETER :: irain = 6  ! indice de traceurs pour gouttes de pluie
106      INTEGER, PARAMETER :: itke  = 7  ! indice de traceurs pour tke
107      INTEGER, PARAMETER :: itked = 8  ! indice de traceurs pour tke dissipation
108! #cw INTEGER, PARAMETER :: iccn  = 9  ! indice de traceurs pour CCN
109
110!======================================================================
111! Arguments de la routine physiq.F90:
112!======================================================================
113      integer,intent(in) :: nlon ! number of atmospheric colums
114      integer,intent(in) :: nlev ! number of vertical levels (should be =klev)
115      real,intent(in) :: jD_cur ! current day number (Julian day)
116      real,intent(in) :: jH_cur ! current time of day (as fraction of day)
117      logical,intent(in) :: debut ! signals first call to physics
118      logical,intent(in) :: lafin ! signals last call to physics
119      real,intent(in) :: pdtphys ! physics time step (s)
120      real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa)
121      real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa)
122      real,intent(in) :: pphi(klon,klev) ! geopotential at mid-layer
123      real,intent(in) :: pphis(klon) ! surface geopotential
124      real,intent(in) :: ppresnivs(klev) ! pseudo-pressure (Pa) of mid-layers
125      real,intent(in) :: u(klon,klev) ! eastward zonal wind (m/s)
126      real,intent(in) :: v(klon,klev) ! northward meridional wind (m/s)
127      REAL, intent(in):: rot(klon, klev) ! relative vorticity,
128                                         ! in s-1, needed for frontal waves
129      real,intent(in) :: t(klon,klev) ! temperature (K)
130      real,intent(in) :: qx(klon,klev,nqtot) ! tracers (.../kg_air)
131      real,intent(in) :: flxmass_w(klon,klev) ! vertical mass flux
132      real,intent(out) :: d_u(klon,klev) ! physics tendency on u (m/s/s)
133      real,intent(out) :: d_v(klon,klev) ! physics tendency on v (m/s/s)
134      real,intent(out) :: d_t(klon,klev) ! physics tendency on t (K/s)
135      real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers
136      real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure
137      real,intent(in) :: dudyn(iim+1,jjmp1,klev) ! Not used
138
139!======================================================================
140! Variables locales temporelles
141!======================================================================
142integer,save :: itau=0 ! counter to count number of calls to physics
143!$OMP THREADPRIVATE(itau)
144!real :: temp_newton(klon,klev)
145integer :: k
146logical, save :: first=.true.
147!$OMP THREADPRIVATE(first)
148
149!======================================================================
150! Inputs/Outputs
151!======================================================================
152integer :: itau0
153real :: zjulian
154real :: dtime
155integer :: nhori ! horizontal coordinate ID
156integer,save :: nid_hist ! output file ID
157!$OMP THREADPRIVATE(nid_hist)
158integer :: zvertid ! vertical coordinate ID
159integer,save :: iwrite_phys ! output every iwrite_phys physics step
160!$OMP THREADPRIVATE(iwrite_phys)
161integer :: iwrite_phys_omp ! intermediate variable to read iwrite_phys
162real :: t_ops ! frequency of the IOIPSL operations (eg average over...)
163real :: t_wrt ! frequency of the IOIPSL outputs
164
165
166!======================================================================
167! Arguments d'entrée pour la physique de MAR
168!======================================================================
169
170! NB : Pour plus de visibilité pour les développeurs de LMDZ et MAR, on conserve les noms de variables propres à MAR et LMDZ.
171! La correspondance entre les variables se fait dans cette interface.
172
173! Flags pour MAR (A ranger ultérieurement dans un fichier .def, selon le protocole LMDZ).
174! Attention, ces Flags pilotent l'initialisation de nombreuses variables. Toutes les configurations n'ont pas été testées.
175! Ici, on trouve les déclarations utilisées pour une experience aquaplanète.
176! Contacter Hubert Gallée pour connaître la signification des Flags.
177
178
179       logical                    ::  FlagSV     = .TRUE.               ! Surface
180       logical                    ::  FlagSV_Veg = .TRUE.               ! Végétation
181       logical                    ::  FlagSV_SNo = .TRUE.               ! Neige
182       logical                    ::  FlagSV_BSn = .FALSE.              !
183       logical                    ::  FlagSV_KzT = .TRUE.               !
184       logical                    ::  FlagSV_SWD = .TRUE.               !  (T,F) == SW INPUT is (downward,absorbed)
185       ! FlagSV_SBC est utilisé pour lire des options de conditions aux limites. Ici, on le met à FALSE pour l'aquaplanète
186       logical                    ::  FlagSV_SBC = .FALSE.              ! Flag pour utiliser le PHY_SISVAT_INP propre à notre experience
187       logical                    ::  FlagSV_UBC = .FALSE.              !
188       logical                    ::  FlagAT     = .TRUE.               !
189       character(len=1)           ::  TypeAT     = 'e'                  !  Turbulence Model (e, K, L or H)
190       logical                    ::  FlagAT_TKE = .TRUE.               !
191       logical                    ::  FlagCM     = .TRUE.               !
192!Gilles : T > F car LMDz utilise les tendances des traceurs (sauf frac.nuage)
193       logical                    ::  FlagCM_UpD = .FALSE.              ! Update de la microphysique immédiat
194       logical                    ::  FlagCP     = .TRUE.               ! Avec ou sans convection
195       logical                    ::  FlagRT     = .TRUE.               !
196       logical                    ::  FlagS0_SLO = .FALSE.              !
197       logical                    ::  FlagS0_MtM = .FALSE.              !
198       logical                    ::  Flag_O     = .TRUE.               !
199!Gilles: FlagVR  T > F pas de sortie interessante
200       logical                    ::  FlagVR     = .FALSE.              !
201
202!======================================================================
203! Arguments pour la physique de MAR :
204!======================================================================
205
206! Pas de temps par défaut. Ils sont modifiés plus bas en fonction de pdtphys défini dans les fichiers .def.
207
208       real                       ::  dt0DYn     =   60.                !  Time Step, Dynamics                       [s] recalculé en fonction de day_step
209       real                       ::  dt0_SV     =   60.                !  Time Step, SISVAT                         [s] f(pdtphys)
210       real                       ::  dt0_AT     =   60.                !  Time Step, Atm_AT                         [s] f(pdtphys)
211       real                       ::  dt0_CM     =   60.                !  Time Step, CMiPhy                         [s] f(pdtphys)
212       real                       ::  dt0_CP     =  600.                !  Time Step, CVAmnh                         [s] f(pdtphys)
213       real                       ::  dt0_RT     =  900.                !  Time Step, radCEP                         [s] f(pdtphys)
214       real                       ::  dx         = 4.0e+4               !  Grid Spacing, MAX authorised for CMiPhy   [m] Limite max tolérée pour la microphysique MAR (40)
215       real(kind=real8), dimension(klon)      ::  dx___HOST             !  Grid Spacing max along x axis.
216       real(kind=real8), dimension(klon)      ::  dy___HOST             !  Grid   Size                               [m]
217       real                       ::  DD_AxX =   90.00                  !  x-Axis Direction                          [degree] Standard (90dg = Eastward) [dg] utilise juste pour les conditions aux limite et la position du soleil.
218       real, dimension(klevp1)    ::  s_HOST                            ! MAR ne fonctionne que avec des coordonnées sigma: sigma=(p-Psommet)/(psurface -Psommet)
219       real, dimension(klon)      ::  sh___HOST                         !  Surface Height                            [m] =0 en aquaplanète. A calculer ultérieurement en fonction de l'environnement LMDZ.
220       real(kind=real8), dimension(klon)         ::  sh_a_HOST          !  Topography Anomaly                                        [m]
221       real(kind=real8), dimension(klon)         ::  slopxHOST          !  Slope, x-direction                                        [-]
222       real(kind=real8), dimension(klon)         ::  slopyHOST          !  Slope, y-direction                                        [-]
223       real(kind=real8), dimension(klon)         ::  slopeHOST          !  Slope                                                     [-]
224       real(kind=real8), allocatable             ::  MMaskHOST(:,:)     !  Mountain Mask                                             [-]
225       real, dimension(klon)                     ::  lonh_HOST          !  Longitude                                              [hour]
226       real, dimension(klon)                     ::  latr_HOST          !  Latitude                                             [radian]
227       real                                      ::  ptop_HOST          !  Pressure Model Top                                      [kPa] 0 in LMDZ
228       real, dimension(klon,klevp1)              ::  pkta_HOST          !  Reduced  Potential Temperature                         [KX/s] !  Temperature * p ** (R/Cp) [KX]
229       real, dimension(klon)                     ::  psa__HOST          !  Pressure Thickness                                      [kPa] a recalculer en fonction du geopotentiel de surface: P*=Psurface-Ptop
230       real, dimension(klon,klevp1)              ::  gZa__HOST          !  Geopotential Height                                   [m2/s2]  pphi dans LMDZ
231       real, dimension(klon,klevp1)              ::  gZam_HOST          !  Geopotential Height, mid-level                        [m2/s2] On fait la moyenne entre couche du dessus et couche du dessous.
232       real, dimension(klon,klev)                ::  Ua___HOST          !  Wind Speed, x-direction                                 [m/s]
233       real, dimension(klon,klev)                ::  Va___HOST          !  Wind Speed, y-direction                                 [m/s]
234       real, dimension(klon,klev)                ::  Wa___HOST          !  Wind Speed, z-direction                                 [m/s] Non echangé entre phys et dyn dans LMDZ. Necessaire pour MAR. On le recalcule ci-aprés !
235       real, dimension(klon,klevp1)              ::  qv___HOST          !  Specific Humidity                                     [kg/kg] qx(:,:,ivap)
236       real, dimension(klon,klev)                ::  qw___HOST          !  Cloud Droplets Concentration                          [kg/kg] qx(:,:,iliq)
237! #cw  real, dimension(klon,klev)                ::  CCN__HOST          !  CCN            Concentration                           [-/kg] qx(:,:,iccn)
238       real, dimension(klon,klev)                ::  qi___HOST          !  Cloud Crystals Concentration                          [kg/kg] qx(:,:,iice)
239       real, dimension(klon,klev)                ::  CIN__HOST          !  CIN            Concentration                           [-/kg] qx(:,:,icin)
240       real, dimension(klon,klev)                ::  CF___HOST          !  Cloud Fraction                                         [-/kg] Dans le futur, la CF sera transportée dans la dynamique, mais pour l'instant, on n'en a pas besoin.
241
242       real, dimension(klon,klev)                ::  qs___HOST          !  Snow Particles Concentration                          [kg/kg] qx(:,:,isnow)
243       real, dimension(klon,klev)                ::  qr___HOST          !  Rain Drops     Concentration                          [kg/kg] qx(:,:,irain)
244       real, dimension(klon,klev)                ::  TKE__HOST          !  Turbulent Kinetic Energy                              [m2/s2] qx(:,:,itke)
245       real, dimension(klon,klev)                ::  eps__HOST          !  Turbulent Kinetic Energy Dissipation                  [m2/s3] qx(:,:,itked)
246       real, dimension(klon,klev)                ::  dpkt___dt          !  Reduced  Potential Temperature TENDENCY, ALL Contribut.[KX/s] a recalculer en fonction de d_t de LMDZ
247       real, dimension(klon,klev)                ::  dua____dt          !  Wind Speed       (x-direc.)    TENDENCY, ALL Contribut.[m/s2] output commune MAR-LMD
248       real, dimension(klon,klev)                ::  dva____dt          !  Wind Speed       (y-direc.)    TENDENCY, ALL Contribut.[m/s2] output commune MAR-LMD
249       real, dimension(klon,klev)                ::  dqv____dt          !  Specific           Humidity    TENDENCY, ALL Contr. [kg/kg/s] output commune MAR-LMD
250       real, dimension(klon,klev)                ::  dqw____dt          !  Cloud Droplets Concentration   TENDENCY, ALL Contr. [kg/kg/s] dqx
251! #cw real, dimension(klon,klev)                 ::  dCw____dt          !  CCN            Concentration   TENDENCY, ALL Contr.  [1/kg/s] dqx
252       real, dimension(klon,klev)                ::  dqi____dt          !  Cloud Crystals Concentration   TENDENCY, ALL Contr. [kg/kg/s] dqx
253       real, dimension(klon,klev)                ::  dCi____dt          !  CIN            Concentration   TENDENCY, ALL Contr.  [1/kg/s] dqx
254       real, dimension(klon,klev)                ::  dCF____dt          !  Cloud Fraction                 TENDENCY, ALL Contr. [kg/kg/s] dqx
255       real, dimension(klon,klev)                ::  dqs____dt          !  Snow Particles Concentration   TENDENCY, ALL Contr. [kg/kg/s] dqx
256       real, dimension(klon,klev)                ::  dqr____dt          !  Rain Drops     Concentration   TENDENCY, ALL Contr. [kg/kg/s] dqx
257! #dT& real, dimension(klon,klev)                ::  dpktSV_dt          !  Reduced  Potential Temperature Tendency, SISVAT        [KX/s]
258! #dT       real, dimension(klon,klev)           ::  dpktAT_dt          !  Reduced  Potential Temperature Tendency, Atm_AT        [KX/s]
259! #dT       real, dimension(klon,klev)           ::  dqv_AT_dt          !  Specific           Humidity    TENDENCY, Atm_AT     [kg/kg/s] dqx
260! #dT       real, dimension(klon,klev)           ::  dqw_AT_dt          !  Cloud Droplets Concentration   TENDENCY, Atm_AT     [kg/kg/s] dqx
261! #cw real, dimension(klon,klev)                 ::  dCw_AT_dt          !  CCN            Concentration   TENDENCY, Atm_AT         [1/s] dqx
262! #dT       real, dimension(klon,klev)           ::  dqi_AT_dt          !  Cloud Crystals Concentration   TENDENCY, Atm_AT     [kg/kg/s] dqx
263! #dT      real, dimension(klon,klev)            ::  dCi_AT_dt          !  CIN            Concentration   TENDENCY, Atm_AT         [1/s] dqx
264! #dT      real, dimension(klon,klev)            ::  dqs_AT_dt          !  Snow Particles Concentration   TENDENCY, Atm_AT     [kg/kg/s] dqx
265! #dT      real, dimension(klon,klev)            ::  dqr_AT_dt          !  Rain Drops     Concentration   TENDENCY, Atm_AT     [kg/kg/s] dqx
266! #dT      real, dimension(klon,klev)            ::  dpktCM_dt          !  Reduced  Potential Temperature Tendency, CMiPhy        [KX/s] dqx
267! #dT      real, dimension(klon,klev)            ::  dqv_CM_dt          !  Specific           Humidity    TENDENCY, CMiPhy     [kg/kg/s] dqx
268! #dT      real, dimension(klon,klev)            ::  dqw_CM_dt          !  Cloud Droplets Concentration   TENDENCY, CMiPhy     [kg/kg/s] dqx
269! #cw real, dimension(klon,klev)                 ::  dCw_CM_dt          !  CCN            Concentration   TENDENCY, CMiPhy         [1/s] dqx
270! #dT       real, dimension(klon,klev)           ::  dqi_CM_dt          !  Cloud Crystals Concentration   TENDENCY, CMiPhy     [kg/kg/s] dqx
271! #dT       real, dimension(klon,klev)           ::  dCi_CM_dt          !  CIN            Concentration   TENDENCY, CMiPhy         [1/s] dqx
272! #dT       real, dimension(klon,klev)           ::  dCF_CM_dt          !  Cloud Fraction                 TENDENCY, CMiPhy         [1/s] dqx
273! #dT       real, dimension(klon,klev)           ::  dqs_CM_dt          !  Snow Particles Concentration   TENDENCY, CMiPhy     [kg/kg/s] dqx
274! #dT       real, dimension(klon,klev)           ::  dqr_CM_dt          !  Rain Drops     Concentration   TENDENCY, CMiPhy     [kg/kg/s] dqx
275! #dT       real, dimension(klon,klev)           ::  dpktCP_dt          !  Reduced  Potential Temperature TENDENCY, CVAmnh        [KX/s]
276! #dT       real, dimension(klon,klev)           ::  dqv_CP_dt          !  Specific           Humidity    TENDENCY, CVAmnh     [kg/kg/s] dqx
277! #dT       real, dimension(klon,klev)           ::  dqw_CP_dt          !  Cloud Droplets Concentration   TENDENCY, CVAmnh     [kg/kg/s] dqx
278! #dT       real, dimension(klon,klev)           ::  dqi_CP_dt          !  Cloud Crystals Concentration   TENDENCY, CVAmnh     [kg/kg/s] dqx
279! #dT       real, dimension(klon,klev)           ::  dpktRT_dt          !  Reduced  Potential Temperature TENDENCY, radCEP        [KX/s]
280
281       integer,save                              ::  it0EXP             ! Incremente en fonction de debut et lafin. NB: MAR fait des restarts dans le Fortran, LMDZ avec des scripts shells. On suit ici la logique LMDZ. it0EXP et it0RUN sont donc identiques.
282       integer,save                              ::  it0RUN             ! Incremente en fonction de debut et lafin.
283       integer                                   ::  Year_H             !  Time                                                   [year]
284       integer                                   ::  Mon__H             !  Time                                                  [month]
285       integer                                   ::  Day__H             !  Time                                                    [Day]
286       integer                                   ::  Hour_H             !  Time                                                   [hour]
287       integer                                   ::  minu_H             !  Time                                                 [minute]
288       integer                                   ::  sec__H             !  Time                                                      [s]
289       integer                                   ::  ixq1,i0x0,mxqq     !  Domain  Dimension: x                                      [-]
290       integer                                   ::  jyq1,j0y0,myqq     !  Domain  Dimension: y                                      [-]
291       !integer                                   ::  klev              !  Domain  Dimension: z                                      [-] Commun à MAR et LMDZ
292       !integer                                   ::  klevp1            !  Domain  Dimension: z                                      [-] Commun à MAR et LMDZ
293       integer                                   ::  mwq                !  Domain  Dimension: mosaic                                 [-]
294       !integer                                   ::  klon              !  Domain  Dimension: x * y                                  [-] Commun à MAR et LMDZ
295       integer                                   ::  kcolw              !  Domain  Dimension: x * y * mosaic                         [-]
296       integer                                   ::  m_azim             !  Mountain Mask, nb of directions taken into account        [-]
297       integer,parameter                         ::  n0pt=1             !  nombre de points qui peuvent être choisis pour le fichier PHY____.OUT
298       integer, dimension(n0pt)                  ::  IOi0SV             !
299       integer, dimension(n0pt)                  ::  IOj0SV             !
300       real, dimension(klon)                     ::  sst__HOST          !  Ocean     FORCING (SST)                                   [K]
301! #IP real, dimension(klon)                       ::  sif__HOST          !  Ocean     FORCING (Sea-Ice Fraction )                     [-]
302! #AO real, dimension(ixq1:mxqq,jyq1:myqq,mwq)    ::  s_T__HOST          !  A - O    COUPLING                      n=1: Open Ocean    [K]
303! #AO real, dimension(ixq1:mxqq,jyq1:myqq,mwq)    ::  Alb__HOST          !  A - O    COUPLING (Surface Albedo   )  n=2: Sea  Ice      [-]
304! #AO real, dimension(ixq1:mxqq,jyq1:myqq,mwq)    ::  dSdT2HOST          !  A - O    COUPLING ( d(SH Flux) / dT )                [W/m2/K]
305! #AO real, dimension(ixq1:mxqq,jyq1:myqq,mwq)    ::  dLdT2HOST          !  A - O    COUPLING ( d(SH Flux) / dT )                [W/m2/K]
306! #TC integer, parameter                          ::   ntrac  =  28      !
307! #TC real, dimension(ixq1:mxqq,jyq1:myqq,klev,ntrac) ::    qxTC             !  Aerosols: Atmospheric  Contentration
308! #TC real, dimension(ixq1:mxqq,jyq1:myqq    ,ntrac) ::    qsTC             !  Aerosols: Near Surface Contentration
309! #TC real, dimension(ixq1:mxqq,jyq1:myqq    ,ntrac) ::    uqTC             !  Aerosols: Surf.Flux
310
311!======================================================================
312! Variables locales:
313!======================================================================
314
315        real :: secondes ! pour le calcul de la date
316        integer :: i ! Pour les boucles
317        integer :: daysec ! Nombre de secondes dans ue journée
318        integer,save :: annee_ref ! Année du début de la simulation
319        !integer :: count_year ! le nombre d'années depuis le début de la simulation
320        !real,parameter :: epsilon = 1.e-6 ! pour le calcul du temps
321
322! Constantes définies comme dans MAR et utilisées ici:
323        real    ::  RdA        =  287.                 !  Perfect Gas Law Constant for Dry Air [J/kg/K]
324        real    ::  CpA        = 1004.                 !  Specific   Heat Capacity for Dry Air [J/kg/K]
325        real    ::  kap                                !  RdA/CpA                                   [-]
326
327! Pour l'initialisation:
328        real, dimension(klon) :: Ts___HOST ! Temperature de surface pour l'initialisation.
329        real, dimension(klon) :: pkt_DY_surf_tmp ! pkt de surface pour l'initialisation.
330        real, dimension(klon) :: qv_HOST_surf_tmp ! qv de surface pour l'initialisation.
331
332! Pour ecriture en netcdf
333       real, dimension(klon,klev) :: cldfra   ! CF___HOST vertically inverted
334
335!======================================================================
336!! Ecriture des sorties netcdf:
337!======================================================================
338
339! NB : La version de MAR utilisée ici ne permet pas d'écrire de sorties netcdf.
340! On n'utilise pas encore XIOS
341! En attendant, on utilise ici le classique ioipsl, ainsi que l'outil iotd de Frédéric Hourdin.
342
343if (debut) then ! Things to do only for the first call to physics
344! load initial conditions for physics (including the grid)
345  call phys_state_var_init() ! some initializations, required before calling phyetat0
346  call phyetat0("startphy.nc")
347
348! Initialize outputs:
349  itau0=0
350!$OMP MASTER
351  iwrite_phys_omp=1 !default: output every physics timestep
352  ! NB: getin() is not threadsafe; only one thread should call it.
353  call getin("iwrite_phys",iwrite_phys_omp)
354!$OMP END MASTER
355!$OMP BARRIER
356  iwrite_phys=iwrite_phys_omp
357  t_ops=pdtphys*iwrite_phys ! frequency of the IOIPSL operation
358  t_wrt=pdtphys*iwrite_phys ! frequency of the outputs in the file
359  ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0
360  !CALL ymds2ju(annee0, month, dayref, hour, zjulian)
361  call getin('anneeref',anneeref)
362  call ymds2ju(anneeref, 1, 1, 0.0, zjulian)
363  annee_ref=anneeref ! On la sauve pour le reste de la simulation
364  dtime=pdtphys
365
366! Initialisation de iophy pour l'utilisation de iotd créé par Frédéric Hourdin
367  call iophys_ini ! NB: cette initialisation concerne l'ecriture de fichiers phys.nc dans PHY_MAR (rechercher iophys_ecrit)
368
369!Gilles : replace zjulian by jD_cur (if correct init of the year)
370! call histbeg_phy("histins.nc",itau0,zjulian,dtime,nhori,nid_hist)
371  call histbeg_phy("histins.nc",itau0,jD_cur,dtime,nhori,nid_hist)
372
373!$OMP MASTER
374
375  ! define vertical coordinate
376  call histvert(nid_hist,"presnivs","Vertical levels","Pa",klev, &
377                ppresnivs,zvertid,'down')
378  ! define variables which will be written in "histins.nc" file
379  call histdef(nid_hist,'temperature','Atmospheric temperature','K', &
380               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
381               'inst(X)',t_ops,t_wrt)
382  call histdef(nid_hist,'pplay','atmospheric pressure','K', &
383               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
384               'inst(X)',t_ops,t_wrt)
385  call histdef(nid_hist,'ps','Surface Pressure','Pa', &
386               iim,jj_nb,nhori,1,1,1,zvertid,32, &
387               'inst(X)',t_ops,t_wrt)
388  call histdef(nid_hist,'qx_vap','specific humidity','kg/kg', &
389               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
390               'inst(X)',t_ops,t_wrt)
391  call histdef(nid_hist,'qx_liq','atmospheric liquid water content','kg/kg', &
392               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
393               'inst(X)',t_ops,t_wrt)
394  call histdef(nid_hist,'qx_ice','atmospheric ice content','kg/kg', &
395               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
396               'inst(X)',t_ops,t_wrt)
397  call histdef(nid_hist,'u','Eastward Zonal Wind','m/s', &
398               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
399               'inst(X)',t_ops,t_wrt)
400  call histdef(nid_hist,'v','Northward Meridional Wind','m/s', &
401               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
402               'inst(X)',t_ops,t_wrt)
403  call histdef(nid_hist,'lonh','longitude','Hours', &
404               iim,jj_nb,nhori,1,1,1,zvertid,32, &
405               'inst(X)',t_ops,t_wrt)
406  call histdef(nid_hist,'latr','latitude','radian', &
407               iim,jj_nb,nhori,1,1,1,zvertid,32, &
408               'inst(X)',t_ops,t_wrt)
409  call histdef(nid_hist,'sst','Prescribed SST','K', &
410               iim,jj_nb,nhori,1,1,1,zvertid,32, &
411               'inst(X)',t_ops,t_wrt)
412  ! Tendances physiques de MAR ; pour le developpement MAR-LMDZ:
413  call histdef(nid_hist,'d_t','Atmospheric temperature trends after MAR physics','K/s', &
414               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
415               'inst(X)',t_ops,t_wrt)
416  call histdef(nid_hist,'d_u','Eastward Zonal Wind trend','m/s/s', &
417               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
418               'inst(X)',t_ops,t_wrt)
419  call histdef(nid_hist,'d_v','Northward Meridional Wind trend','m/s/s', &
420               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
421               'inst(X)',t_ops,t_wrt)
422  ! Définition de variables utilisée dans des sous-routines MAR
423  call histdef(nid_hist,'snow','snowfall convectif + stratiform','m', &
424               iim,jj_nb,nhori,1,1,1,zvertid,32, &
425               'inst(X)',t_ops,t_wrt)
426  call histdef(nid_hist,'rain','rainfall convectif + stratiform','m', &
427               iim,jj_nb,nhori,1,1,1,zvertid,32, &
428               'inst(X)',t_ops,t_wrt)
429  call histdef(nid_hist,'snowCP','snowfall convective','m', &
430               iim,jj_nb,nhori,1,1,1,zvertid,32, &
431               'inst(X)',t_ops,t_wrt)
432  call histdef(nid_hist,'rainCP','rainfall convective','m', &
433               iim,jj_nb,nhori,1,1,1,zvertid,32, &
434               'inst(X)',t_ops,t_wrt)
435  call histdef(nid_hist,'nebulosity','CF___HOST','-', &
436               iim,jj_nb,nhori,klev,1,klev,zvertid,32, &
437               'inst(X)',t_ops,t_wrt)
438  ! end definition sequence
439  call histend(nid_hist)
440
441!XIOS
442#ifdef CPP_XIOS
443    ! Déclaration de l'axe vertical du fichier:   
444    !CALL wxios_add_vaxis("presnivs", "histins", klev, ppresnivs)
445
446    !Déclaration du pas de temps:
447    !CALL wxios_set_timestep(dtime)
448
449    !Finalisation du contexte:
450    !CALL wxios_closedef()
451#endif
452!$OMP END MASTER
453endif ! of if (debut)
454
455!======================================================================
456! Calcul de variables environnement nécessaires à la physique de MAR
457!======================================================================
458
459! NB : Calculs à faire avant l'initialisation de la physique
460
461! Calcul de constantes
462
463kap    = RdA  / CpA
464rpi  = 4.*ATAN(1.)
465daysec = 86400.
466
467! increment counter itau
468itau=itau+1
469
470PRINT*, 'dans physiq.F90 :'
471PRINT*, 'itau=',itau
472
473! set all tendencies to zero
474PRINT*, 'On initialise les tendances à zéro:'
475d_u(1:klon,1:klev)=0.
476d_v(1:klon,1:klev)=0.
477d_t(1:klon,1:klev)=0.
478d_qx(1:klon,1:klev,1:nqtot)=0.
479d_ps(1:klon)=0.
480
481! Allocation de tableaux:
482! Pour l'instant :
483m_azim=1 ! Azimuths pour le calcul des masques de montagne
484ALLOCATE(MMaskHOST(klon,m_azim)) !  Mountain Mask                                             [-]
485
486! Appeller iniaqua et phyetat0 pour l'état initial avant le transfert vers les variables MAR !
487! A faire pour lorsque l'on prendra en compte les continents
488
489! Timesteps pour MAR (définis ici à partir des pas de temps LMDZ). A modifier éventuellement (via les fichiers .def)
490! dt0_CP & dt0_RT modifiees ici selon Hubert avec daystep=720 ie dtdyn=2mn
491! pdtphys modifie via gcm.def avec iphysiq=1 & nsplit_phys=1
492dt0DYn=daysec/REAL(day_step)   !  Time Step, Dynamics                       [s]
493dt0_AT=pdtphys                 !  Time Step, SISVAT                         [s]
494dt0_SV=pdtphys                 !  Time Step, Atm_AT                         [s]
495dt0_CM=pdtphys                 !  Time Step, CMiPhy                         [s]
496dt0_CP=pdtphys*10              !  Time Step, CVAmnh                         [s]
497dt0_RT=pdtphys*20              !  Time Step, radCEP                         [s]
498
499! Correspondances de grilles: inutilisé car on utilise directement les variables LMD
500
501! Dans un premier temps, une seule mozaique:
502mwq=1
503
504kcolw=mwq*klon
505
506DO i=1,klon
507  IF (longitude(i) .LT. 0) THEN
508    lonh_HOST(i)=longitude(i)*12./rpi+24. ! from radians to hours
509   ELSE
510    lonh_HOST(i)=longitude(i)*12./rpi ! from radians to hours
511  ENDIF
512ENDDO
513latr_HOST(:)=latitude(:) ! from radians to radians
514
515!PRINT*, 'lonh_HOST(:)=',lonh_HOST(:)
516!PRINT*, 'latr_HOST(:)=',latr_HOST(:)
517
518! Niveaux verticaux:
519! Attention, il faut faire un test pour vérifier si les coordonnées verticaux définis sont adéquats avec ceux de MAR!
520DO k = 1, klev
521  i=klev+1-k
522  IF (ap(k) .NE. 0) THEN
523    CALL abort_gcm('','MAR physic can be applied only with sigma levels, check vertical resolution in disvert.F90',1)
524   ELSE
525    s_HOST(i)=0.5*(bp(k)+bp(k+1))
526  ENDIF
527END DO
528s_HOST(klev+1)=1 ! Dans MAR, sigma est défini depuis tout en haut jusqu'à la surface. s_HOST doit être égal à 1 en surface.
529
530! En attendant que les routines PHY_____INI.F90, PHY_MAR_DAT.F90 et PHY_SISVAT_INP.F90 soient écrites selon un vecteur et non pas selon un tableau 2D lon-lat.
531! On utilise ici la grille physique de LMDZ, qui a la taille klon (en argument de la physique): 2 + (JJM-1)*IIM (soit 1682 pour une grille dynamique de 48*36 déclarée à la compilation).
532
533mxqq=iim
534myqq=jjm+1
535
536! Pour MAR, on a besoin des indices du premier point de la grille. Avec une grille globale, c'esty toujours (1,1):
537ixq1=1
538jyq1=1
539
540
541! Temps:
542IF (debut) THEN
543  it0EXP=1
544 ELSE
545  it0EXP=it0EXP+1
546ENDIF
547
548IF (debut) THEN
549  it0RUN=1
550 ELSE
551  it0RUN=it0RUN+1
552ENDIF
553
554! Pour l'instant, on ne peut démarrer une simulation seulement le 1er jour de l'année de référence:
555! Attention, ju2ymds ne donne pas l'année courante, mais compte les années de simulations à partir de 0 !!! A checker !!!
556CALL ymds2ju(annee_ref, 1, 1, 0.0, zjulian)
557!Gilles: zjulian already counted in jD_cur if year correctly initialised
558!CALL ju2ymds(jD_cur+jH_cur+zjulian,Year_H,Mon__H,Day__H,secondes)
559CALL ju2ymds(jD_cur+jH_cur,Year_H,Mon__H,Day__H,secondes)
560Hour_H=INT(secondes)/3600
561minu_H=(INT(secondes)-Hour_H*3600)/60
562!minu_H=INT((secondes-float(Hour_H)*3600.+epsilon)/60.) ! attention aux arrondis machine
563sec__H=INT(secondes)-(Hour_H*3600)-(minu_H*60)
564
565Print*,'CONTROL du temps avant d appeler la physique'
566Print*,'zjulian=',zjulian
567Print*,'annee_ref=',annee_ref
568Print*, 'jD_cur=',JD_cur
569Print*, 'jH_cur=',jH_cur
570Print*, 'Year_H=',Year_H
571Print*, 'Mon__H',Mon__H
572Print*, 'Day__H',Day__H
573Print*, 'secondes=',secondes
574Print*, 'Hour_H=',Hour_H
575Print*, 'minu_H=',minu_H
576Print*, 'sec__H=',sec__H
577Print*, 'klev=',klev
578Print*, 'klevp1=',klevp1
579Print*, 'nlev=',nlev
580Print*, 'nlon=',nlon
581Print*, 'klon=',klon
582PRINT*, 'RCp=',RCp
583
584! Coordonnées du point (en indice) pour lequel on imprime les sorties en ASCII dans le fichier PHY___________.OUT
585i0x0=23
586j0y0=48
587! Variable utilisée juste pour les outputs:
588IOi0SV(:)=1
589IOj0SV(:)=1
590
591! Variables inutilisees:
592dy___HOST=0.
593!CF___HOST=0. ! Pour l'instant pas utilisé, il faut tout de même l'initialiser?
594dCF____dt=0.
595
596
597!======================================================================
598! INITIALISATIONS MAR
599!======================================================================
600!
601! NB : les routines MAR ci-aprés doivent encore être travaillées. Ici, elles ne fonctionnent que pour une aquaplanète !
602
603       it_RUN = it0RUN
604       it_EXP = it0EXP
605
606       IF (it0RUN.LE.1)                                              THEN
607
608         PRINT*,'it0RUN=0; On appelle les routines d initialisation de la physique de MAR'
609
610! Initialization of Mod_PHY____grd
611! --------------------------------
612
613! Anciennes initialisations MAR
614
615         dxHOST = dx
616         DirAxX = DD_AxX
617         ixp1   = ixq1
618         i_x0   = i0x0
619         mxpp   = mxqq
620         jyp1   = jyq1
621         j_y0   = j0y0
622         mypp   = myqq
623         mzpp   = klevp1
624         kcolp  = klon
625         kcolv  = kcolw
626
627! Initialization of Mod_PHY_S0_grd
628! --------------------------------
629         ! Pour l'instant, en attendant du relief :
630         n_azim = m_azim
631
632! Initialization of Mod_SISVAT_grd (Mosaic Discretization only                                      )
633! -------------------------------- (       Discretization from Mod_SISVAT_dim : see PHY_SISVAT_ALLOC)
634
635          mwp    = mwq
636
637! Initialization of Mod_SISVAT_dat
638! --------------------------------
639! pour aquaplanete: FixedSST=1
640          VarSST = 0.
641
642! Initialization of Physical Constants and Allocation of main Variables
643! ---------------------------------------------------------------------
644
645! CONTROL
646PRINT*, 'Appel de PHY INI'
647!              **************
648         CALL PHY________INI
649!              **************
650
651
652! INPUT DATA surface et Calcul des dependances a des Points Horizontaux autres que le Point courant
653! =================================================================================================
654!
655
656! Pour l'instant, PHY_MAR_DAT n'est pas vectorisé, donc on initilise tout en dur pour l'aquaplanète.
657! Ultérieurement, il faudra lires les données topo de LMDZ pour renseigner ces variables à l'initialisation.
658
659PRINT*,'Initialisation en dur des variables topo pour une aquaplanete'
660sh___HOST(:)=0.00
661sh_a_HOST(:)=0.00
662dx___HOST(:)   = 1.e3
663dy___HOST(:)   = 1.e3
664slopxHOST(:)   = 0.00
665slopyHOST(:)   = 0.00
666slopeHOST(:)   = 0.00
667MMaskHOST(:,:) = 0.00
668
669! CONTROL
670!PRINT*, 'Appel de PHY MAR DAT'
671!          ***********
672!      CALL PHY_MAR_DAT(                                                &
673!!          ***********
674!&                 sh___HOST                                       &!  Topographie                                 [m]
675!&                ,sh_a_HOST                                       &!  Topographie, Anomalie                       [m]
676!&                ,dx___HOST                                       &!  Discretisation,   x                         [m]
677!&                ,dy___HOST                                       &!  Discretisation,   y                         [m]
678!&                ,slopxHOST                                       &!  Pente selon l'Axe x                         [-]
679!&                ,slopyHOST                                       &!  Pente selon l'Axe y                         [-]
680!&                ,slopeHOST                                       &!  Pente                                       [-]
681!&                ,MMaskHOST                                       &!  Masque de Montagnes                         [-]
682!&                ,ixq1,mxqq                                       &!
683!&                ,jyq1,myqq                                       &!
684!&                ,m_azim                                          &!
685!&                )
686
687! FIRST Initialization of the Surface Variables
688! =============================================
689
690
691! Initial Surface Temperature ( initialement dans DYN_to_PHY_MAR)
692! ---------------------------
693
694! CONTROL
695PRINT*, 'Calcul de Ts___HOST dans la physique'
696PRINT*, 'Initialisation de la temperature de surface avec les SST aquaplanète'
697          DO i = 1,klon
698            Ts___HOST(i)=273.+27.*(1-sin(1.5*latitude(i))**2)
699            IF ((latitude(i).GT.1.0471975).OR.(latitude(i).LT.-1.0471975)) THEN
700              Ts___HOST(i)=273.
701            ENDIF
702
703             !PRINT*, 'Initialisation de la temperature de surface avec la temperature du premier niveau'
704             !Ts___HOST(i) = pkta_HOST(i,klevp1)*(psa__HOST(i)+ptop_HOST)**RCp
705             !Ts___HOST(i) = t(i,klev)
706          ENDDO
707
708! CONTROL
709PRINT*, 'Avant PHY SISVAT INP'
710PRINT*,'MAXVAL(Ts___HOST(:)=)',maxval(Ts___HOST(:))
711PRINT*,'MINVAL(Ts___HOST(:)=)',minval(Ts___HOST(:))
712PRINT*, 'Appel de PHY SISVAT INP'
713!              **************
714          CALL PHY_SISVAT_INP(FlagSV_SBC,Ts___HOST,klon)
715!              **************
716
717PRINT*, 'Apres PHY SISVAT INP'
718PRINT*,'MAXVAL(Ts___HOST(:)=)',maxval(Ts___HOST(:))
719PRINT*,'MINVAL(Ts___HOST(:)=)',minval(Ts___HOST(:))
720
721       END IF ! (it0RUN.LE.1)
722
723
724! Horizontal Interactions
725! ======================
726
727
728! Interface: DATA for Stand-Alone Version of the Physical Package
729! ===============================================================
730
731       IF  (.NOT. FlagRT .OR.                                           &
732      &     .NOT. FlagAT .OR.                                           &
733      &     .NOT. FlagCM)                                           THEN
734!     &     .NOT. FlagCP)                                            THEN
735! CONTROL
736PRINT*, 'Appel de PHY SISVAT UBC'
737!                **************
738           CALL   PHY_SISVAT_UBC(FlagRT,FlagAT,FlagCM,FlagCP)
739!                **************
740
741       END IF
742
743PRINT*, 'Apres PHY_SISVAT_UBC'
744Print*, 'mzp=',mzp
745Print*, 'mzpp=',mzpp
746
747!======================================================================
748! Calcul des inputs de MAR à partir de LMDZ
749!======================================================================
750
751! NB : Calculs à faire aprés l'initialisation de la physique du MAR
752
753PRINT*, 'Calcul des inputs de MAR à partir de LMDZ'
754
755! Correspondance ou transformation de variables:
756! RKAPPA=R/Cp dans un module LMD ;
757! Probleme, on n'a pas la temperature de surface dans les arguments !!! a voir avec Hubert comment on fait...
758!  En attendant, on prend la température du premier niveau au lieu de celle de la surface (idem pour le haut de l'atmosphère):
759! NB : Attention aux unités !
760
761IF (debut) THEN
762! On initialise avec la temperature du premier niveau au lieu de la temperature de surface à la première itération:
763   pkt_DY_surf_tmp(:)=t(:,1) / (pplay(:,1)/1000) ** kap
764   ELSE
765   pkt_DY_surf_tmp(:)=pkt_DY(:,klev+1)
766ENDIF
767
768pkta_HOST(:,klev+1)=pkt_DY_surf_tmp(:)
769DO k = 1, klev
770  i=klev+1-k
771  pkta_HOST(:,i) = t(:,k) / (pplay(:,k)/1000) ** kap
772END DO
773
774ptop_HOST=0.
775psa__HOST(:)=(paprs(:,1)-paprs(:,klev+1))/1000. ! [kPa]
776
777! Correspondance du géopotentiel entre LMDZ et MAR :
778
779! Aux niveaux:
780DO k = 1, klev
781  i=klev+1-k
782  gZa__HOST(:,i)=pphi(:,k)
783END DO
784gZa__HOST(:,klev+1)=pphis(:)
785
786! Aux inter-niveaux:
787DO k = 1, klev
788  gZam_HOST(:,k+1)=0.5*(gZa__HOST(:,k+1)+gZa__HOST(:,k))
789END DO
790gZam_HOST(:,1)=1.5*gZa__HOST(:,1)-0.5*gZa__HOST(:,2) ! Extrapolation au sommet de l'atmosphère
791DO k = 1, klev
792  i=klev+1-k
793  Ua___HOST(:,i)=u(:,k)
794  Va___HOST(:,i)=v(:,k)
795ENDDO
796
797! Calcul de vitesse verticale a partir de flux de masse verticale :
798DO k = 1, klev
799  i=klev+1-k
800  !omega(i,k) = RG*flxmass_w(i,k) / cell_area(i) ! omega en Pa/s
801  Wa___HOST(:,i)=flxmass_w(:,k) / cell_area(:) * (gZam_HOST(:,i+1) - gZam_HOST(:,i))/(paprs(:,k+1)-paprs(:,k)) ! Equilibre hydrostatique
802END DO
803Wa___HOST(:,klev)=0 ! Vitesse nulle à la surface.
804Wa___HOST(:,1)=0 ! Vitesse nulle au sommet.
805
806! Tendances:
807DO k = 1, klev
808  i=klev+1-k
809  dpkt___dt(:,i)=d_t(:,k) / (pplay(:,k)/1000) ** kap
810  dua____dt(:,i)=d_u(:,k)
811  dva____dt(:,i)=d_v(:,k)
812ENDDO
813
814
815DO k = 1, klev
816  i=klev+1-k
817
818! Traceurs:
819  qv___HOST(:,i)=qx(:,k,ivap)          !  Specific Humidity                                     [kg/kg]
820  qw___HOST(:,i)=qx(:,k,iliq)          !  Cloud Droplets Concentration                          [kg/kg]
821  qi___HOST(:,i)=qx(:,k,iice)          !  Cloud Crystals Concentration                          [kg/kg]
822  CIN__HOST(:,i)=qx(:,k,icin)          !  CIN            Concentration                           [-/kg]
823  qs___HOST(:,i)=qx(:,k,isnow)         !  Snow Particles Concentration                          [kg/kg]
824  qr___HOST(:,i)=qx(:,k,irain)         !  Rain Drops     Concentration                          [kg/kg]
825  TKE__HOST(:,i)=qx(:,k,itke)          !  Turbulent Kinetic Energy                              [m2/s2]
826  eps__HOST(:,i)=qx(:,k,itked)         !  Turbulent Kinetic Energy Dissipation                  [m2/s3]
827  ! #cw CCN__HOST(:,:)=qx(:,:,iccn)    !  CCN            Concentration                           [-/kg]
828
829END DO
830
831! Dans MAR, qv(klev+1)=qv(klev)
832IF (debut) THEN
833  ! On initialise avec le qv du premier niveau au lieu du qv de surface à la première itération:
834    qv_HOST_surf_tmp(:)=qv___HOST(:,klev)
835   ELSE
836    qv_HOST_surf_tmp(:)=qv__DY(:,klev+1)
837ENDIF
838qv___HOST(:,klev+1)=qv_HOST_surf_tmp(:)  !  Specific Humidity                                     [kg/kg]
839
840DO k = 1, klev
841  i=klev+1-k
842
843! Dérivée des traceurs:
844dqv____dt(:,i)=d_qx(:,k,ivap)        !  Specific           Humidity    TENDENCY, ALL Contr. [kg/kg/s]
845dqw____dt(:,i)=d_qx(:,k,iliq)        !  Cloud Droplets Concentration   TENDENCY, ALL Contr. [kg/kg/s]
846dqi____dt(:,i)=d_qx(:,k,iice)        !  Cloud Crystals Concentration   TENDENCY, ALL Contr. [kg/kg/s]
847dCi____dt(:,i)=d_qx(:,k,icin)        !  CIN            Concentration   TENDENCY, ALL Contr.  [1/kg/s]
848dqs____dt(:,i)=d_qx(:,k,isnow)       !  Snow Particles Concentration   TENDENCY, ALL Contr. [kg/kg/s]
849dqr____dt(:,i)=d_qx(:,k,irain)       !  Rain Drops     Concentration   TENDENCY, ALL Contr. [kg/kg/s]
850! #cw dCw____dt                      !  CCN            Concentration   TENDENCY, ALL Contr.  [1/kg/s]
851
852END DO
853
854!======================================================================
855! Test Aquaplanète:
856!======================================================================
857
858! CONTROL
859PRINT*, 'Calcul des SST en aquaplanète pour MAR comme dans le cas n°101 de LMDZ' ! (cf iniaqua de LMDZ)
860!IF (debut) THEN
861  DO i=1,klon
862    sst__HOST(i)=273.+27.*(1-sin(1.5*latitude(i))**2)
863    IF ((latitude(i).GT.1.0471975).OR.(latitude(i).LT.-1.0471975)) THEN
864      sst__HOST(i)=273.
865    ENDIF
866  ENDDO
867!ENDIF
868
869!======================================================================
870! CONTROL
871!======================================================================
872PRINT*, 'Control variables LMDZ entree de la physique'
873PRINT*,'MAXVAL(pplay(:,klev))=',MAXVAL(pplay(:,klev))
874PRINT*,'MINVAL(pplay(:,klev))=',MINVAL(pplay(:,klev))
875PRINT*,'MAXVAL(t(:,klev))=',MAXVAL(t(:,klev))
876PRINT*,'MINVAL(t(:,klev))=',MINVAL(t(:,klev))
877PRINT*,'MAXVAL(pplay(:,1))=',MAXVAL(pplay(:,1))
878PRINT*,'MINVAL(pplay(:,1))=',MINVAL(pplay(:,1))
879PRINT*,'MAXVAL(t(:,1))=',MAXVAL(t(:,1))
880PRINT*,'MINVAL(t(:,1))=',MINVAL(t(:,1))
881PRINT*,'MAXVAL(pplay)=',MAXVAL(pplay)
882PRINT*,'MINVAL(pplay)=',MINVAL(pplay)
883PRINT*,'MAXVAL(t)=',MAXVAL(t)
884PRINT*,'MINVAL(t)=',MINVAL(t)
885PRINT*,'MAXVAL(pphis)=',MAXVAL(pphis)
886PRINT*,'MINVAL(pphis)=',MINVAL(pphis)
887PRINT*,'kap=',kap
888PRINT*, 'CONTROL variables MAR entree de la physique'
889PRINT*,'MAXVAL(pkta_HOST)=',MAXVAL(pkta_HOST)
890PRINT*,'MAXVAL(psa__HOST)=',MAXVAL(psa__HOST)
891PRINT*,'MAXVAL(gZa__HOST)=',MAXVAL(gZa__HOST)
892PRINT*,'MAXVAL(gZam_HOST)=',MAXVAL(gZam_HOST)
893PRINT*,'MAXVAL(Ua___HOST)=',MAXVAL(Ua___HOST)
894PRINT*,'MAXVAL(Va___HOST)=',MAXVAL(Va___HOST)
895PRINT*,'MAXVAL(qv___HOST)=',MAXVAL(qv___HOST)
896PRINT*,'MINVAL(pkta_HOST)=',MINVAL(pkta_HOST)
897PRINT*,'MINVAL(psa__HOST)=',MINVAL(psa__HOST)
898PRINT*,'MINVAL(gZa__HOST)=',MINVAL(gZa__HOST)
899PRINT*,'MINVAL(gZam_HOST)=',MINVAL(gZam_HOST)
900PRINT*,'MINVAL(Ua___HOST)=',MINVAL(Ua___HOST)
901PRINT*,'MINVAL(Va___HOST)=',MINVAL(Va___HOST)
902PRINT*,'MINVAL(qv___HOST)=',MINVAL(qv___HOST)
903PRINT*,'Control iterations:'
904PRINT*,'it0RUN=',it0RUN
905PRINT*,'it0EXP=',it0EXP
906PRINT*,'it_RUN=',it_RUN
907PRINT*,'it_EXP=',it_EXP
908
909! On choisi l'endroit où on fait les diagnostics de verification MAR (fichier ASCII):
910! Pour la phase de débug, on cherche des Nan !
911
912DO i=1,klon
913  DO k=1,klev
914    IF (isnan(Va___HOST(i,k))) THEN
915      ikl0 = i
916      PRINT*,'Attention : NaN at'
917      PRINT*,'longitude=',longitude(i)*180/rpi
918      PRINT*,'latitude=',latitude(i)*180/rpi
919    ENDIF
920  ENDDO
921ENDDO
922
923!======================================================================
924! Appel de PHY_MAR
925!======================================================================
926
927! NB: Logiquement, il faut conserver cet appel tel quel, car Hubert Gallée s'engage à ne pas le changer lors des mises à jour !
928
929! CONTROL
930PRINT*, 'Appel de PHY MAR'
931!PRINT*, 'FlagSV     =',FlagSV
932!PRINT*, 'FlagSV_Veg =',FlagSV_Veg
933!PRINT*, 'FlagSV_SNo =',FlagSV_SNo
934!PRINT*, 'FlagSV_BSn =',FlagSV_BSn
935!PRINT*, 'FlagSV_KzT =',FlagSV_KzT
936!PRINT*, 'FlagSV_SWD =',FlagSV_SWD
937!PRINT*, 'FlagSV_SBC =',FlagSV_SBC
938!PRINT*, 'FlagSV_UBC =',FlagSV_UBC
939!PRINT*, 'FlagAT     =',FlagAT
940
941      CALL  PHY_MAR(FlagSV                                        &   ! FLAG  for SISVAT: (T,F) =                    (active OR NOT)
942&                  ,FlagSV_Veg                                    &   ! FLAG  for SISVAT: (T,F) =       (Variable Vegetation OR NOT)
943&                  ,FlagSV_SNo                                    &   ! FLAG  for SISVAT: (T,F) =         (Snow Model active OR NOT)
944&                  ,FlagSV_BSn                                    &   ! FLAG  for SISVAT: (T,F) = (Blowing Snow Model active OR NOT)
945&                  ,FlagSV_KzT                                    &   ! FLAG  for SISVAT: (T,F) = (pkt Turb.Transfert active OR NOT in SISVAT)
946&                  ,FlagSV_SWD                                    &   ! FLAG  for SISVAT: (T,F) = (Modify SW INPUT->downward OR NOT)     
947&                  ,FlagSV_SBC                                    &   ! FLAG  for SISVAT: (T,F) = (INPUT of Soil & Vege DATA OR NOT in SISVAT)
948&                  ,FlagSV_UBC                                    &   ! FLAG  for SISVAT: (T,F) = (pkt UpperBC is Von Neuman OR NOT in SISVAT)
949&                  ,FlagAT                                        &   ! FLAG  for Atm_AT: (T,F) = (Turbulent Transfer active OR NOT)
950&                  ,TypeAT                                        &   ! TYPE  of  Atm_AT: (e= Ee Duynkerke, K= Ee Kitada, L= EL, H= Ee Huan-R)
951&                  ,FlagAT_TKE                                    &   ! FLAG  for genTKE: (T,F) = (TKE-e     Model    active OR NOT)
952&                  ,FlagCM                                        &   ! FLAG  for CMiPhy: (T,F) = (Cloud Microphysics active OR NOT)
953&                  ,FlagCM_UpD                                    &   ! FLAG  for CMiPhy: (T,F) = (qv & hydrometeors updated OR NOT IN CMiPhy)
954&                  ,FlagCP                                        &   ! FLAG  for Convection Paramet.
955&                  ,FlagRT                                        &   ! FLAG  for Radiative Transfer
956&                  ,FlagS0_SLO                                    &   ! FLAG  for Insolation, Surfac.Slope                 Impact  included  NEW
957&                  ,FlagS0_MtM                                    &   ! FLAG  for Insolation, Surfac.Slope & Mountain Mask Impacts included  NEW
958&                  ,Flag_O                                        &   ! FLAG  for OUTPUT
959&                  ,FlagVR                                        &   ! FLAG  for OUTPUT for VERIFICATION
960&                  ,dt0DYn                                        &   ! Time STEP between 2 CALLs of PHY_MAR                     [s] I, fix
961&                  ,dt0_SV                                        &   ! Time STEP between 2 CALLs of SISVAT                      [s] I, fix
962&                  ,dt0_AT                                        &   ! Time STEP between 2 CALLs of Atm_AT                      [s] I, fix
963&                  ,dt0_CM                                        &   ! Time STEP between 2 CALLs of CMiPhy                      [s] I, fix
964&                  ,dt0_CP                                        &   ! Time STEP between 2 CALLs of CVamnh                      [s] I, fix
965&                  ,dt0_RT                                        &   ! Time STEP between 2 CALLs of radCEP                      [s] I, fix
966&                  ,dx                                            &   ! Grid  Mesh size (Horizontal)                             [m] I, fix
967&                  ,DD_AxX                                        &   ! Grid  x-Axis Direction                              [degree] I, fix
968&                  ,s_HOST                                        &   ! Grid (Vertical)   of HOST (NORMALIZED PRESSURE assumed)  [-] I, fix
969&                  ,sh___HOST                                     &   ! Topography                                               [m] I, fix
970&                  ,sh_a_HOST                                     &   ! Topography Anomaly                                       [m] I, fix  NEW
971&                  ,slopxHOST                                     &   ! Slope, x-direction                                       [-] I, fix  NEW
972&                  ,slopyHOST                                     &   ! Slope, y-direction                                       [-] I, fix  NEW
973&                  ,slopeHOST                                     &   ! Slope                                                    [-] I, fix  NEW
974&                  ,MMaskHOST                                     &   ! Mountain Mask                                            [-] I, fix  NEW
975&                  ,lonh_HOST                                     &   ! Longitude                                             [hour] I, fix
976&                  ,latr_HOST                                     &   ! Latitude                                            [radian] I, fix
977&                  ,pkta_HOST                                     &   ! Reduced Potential Temperature                           [XK] I, O
978&                  ,ptop_HOST                                     &   ! Pressure, Model Top                                    [kPa] I, fix
979&                  ,psa__HOST                                     &   ! Pressure  Thickness                                    [kPa] I
980&                  ,gZa__HOST                                     &   ! Geopotential Height                                  [m2/s2] I
981&                  ,gZam_HOST                                     &   ! Geopotential Height, mid-level                       [m2/s2] I
982&                  ,Ua___HOST                                     &   ! Wind ,  x-Direction                                    [m/s] I
983&                  ,Va___HOST                                     &   ! Wind ,  y-Direction                                    [m/s] I
984&                  ,Wa___HOST                                     &   ! Wind ,  z-Direction                                    [m/s] I
985&                  ,qv___HOST                                     &   ! Specific  Humidity                                   [kg/kg] I, O
986&                  ,qw___HOST                                     &   ! Cloud Droplets Concentration                         [kg/kg] I, O
987! #cw&             ,CCN__HOST                                     &   ! CCN            Concentration                          [-/kg]
988&                  ,qi___HOST                                     &   ! Cloud Crystals Concentration                         [kg/kg] I, O
989&                  ,CIN__HOST                                     &   ! CIN            Concentration                          [-/kg] I, O
990&                  ,CF___HOST                                     &   ! Cloud Fraction                                           [-] I, O
991&                  ,qs___HOST                                     &   ! Snow Particles Concentration                         [kg/kg] I, O
992&                  ,qr___HOST                                     &   ! Rain Drops     Concentration                         [kg/kg] I, O
993&                  ,TKE__HOST                                     &   ! Turbulent Kinetic Energy                             [m2/s2] I, O
994&                  ,eps__HOST                                     &   ! Turbulent Kinetic Energy Dissipation                 [m2/s3] I, O
995&                  ,dpkt___dt                                     &   ! Reduced Potential Temperature TENDENCY, ALL Contribut.[KX/s]    O
996&                  ,dua____dt                                     &   ! Wind Speed       (x-direc.)   TENDENCY, ALL Contribut.[m/s2]    O
997&                  ,dva____dt                                     &   ! Wind Speed       (y-direc.)   TENDENCY, ALL Contribut.[m/s2]    O
998&                  ,dqv____dt                                     &   ! Specific          Humidity    TENDENCY, ALL Contr. [kg/kg/s]    O
999&                  ,dqw____dt                                     &   ! Cloud Droplets Concentration  TENDENCY, ALL Contr. [kg/kg/s]    O
1000! #cw&                  ,dCw____dt                                     &   ! CCN            Concentration  TENDENCY, ALL Contr.     [1/s]
1001&                  ,dqi____dt                                     &   ! Cloud Crystals Concentration  TENDENCY, ALL Contr. [kg/kg/s]    O
1002&                  ,dCi____dt                                     &   ! CIN            Concentration  TENDENCY, ALL Contr.     [1/s]    O
1003&                  ,dCF____dt                                     &   ! Cloud Fraction                TENDENCY, ALL Contr.     [1/s]    O
1004&                  ,dqs____dt                                     &   ! Snow Particles Concentration  TENDENCY, ALL Contr. [kg/kg/s]    O
1005&                  ,dqr____dt                                     &   ! Rain Drops     Concentration  TENDENCY, ALL Contr. [kg/kg/s]   O
1006! #dT&                  ,dpktSV_dt                                     &   ! Reduced Potential Temperature TENDENCY, SISVAT        [KX/s]  (O)
1007! #dT&                  ,dpktAT_dt                                     &   ! Reduced Potential Temperature TENDENCY, Atm_AT        [KX/s]  (O)
1008! #dT&                  ,dqv_AT_dt                                     &   ! Specific          Humidity    TENDENCY, Atm_AT     [kg/kg/s]  (O)
1009! #dT&                  ,dqw_AT_dt                                     &   ! Cloud Droplets Concentration  TENDENCY, Atm_AT     [kg/kg/s]  (O)
1010! #dT&                  ,dqi_AT_dt                                     &   ! Cloud Crystals Concentration  TENDENCY, Atm_AT     [kg/kg/s]  (O)
1011! #dT&                  ,dqs_AT_dt                                     &   ! Snow Particles Concentration  TENDENCY, Atm_AT     [kg/kg/s]  (O)
1012! #dT&                  ,dqr_AT_dt                                     &   ! Rain Drops     Concentration  TENDENCY, Atm_AT     [kg/kg/s]  (O)
1013! #cw&                  ,dCw_AT_dt                                     &   ! CCN            Concentration  TENDENCY, Atm_AT         [1/s]  (O)
1014! #dT&                  ,dCi_AT_dt                                     &   ! CIN            Concentration  TENDENCY, Atm_AT         [1/s]  (O)
1015! #dT&                  ,dpktCM_dt                                     &   ! Reduced Potential Temperature TENDENCY, CMiPhy        [KX/s]  (O)
1016! #dT&                  ,dqv_CM_dt                                     &   ! Specific          Humidity    TENDENCY, CMiPhy     [kg/kg/s]  (O)
1017! #dT&                  ,dqw_CM_dt                                     &   ! Cloud Droplets Concentration  TENDENCY, CMiPhy     [kg/kg/s]  (O)
1018! #dT&                  ,dCF_CM_dt                                     &   ! Cloud Fraction                TENDENCY, CMiPhy         [1/s]  (O)
1019! #dT&                  ,dqi_CM_dt                                     &   ! Cloud Crystals Concentration  TENDENCY, CMiPhy     [kg/kg/s]  (O)
1020! #dT&                  ,dqs_CM_dt                                     &   ! Snow Particles Concentration  TENDENCY, CMiPhy     [kg/kg/s]  (O)
1021! #dT&                  ,dqr_CM_dt                                     &   ! Rain Drops     Concentration  TENDENCY, CMiPhy     [kg/kg/s]  (O)
1022! #cw&                  ,dCw_CM_dt                                     &   ! CCN            Concentration  TENDENCY, CMiPhy         [1/s]  (O)
1023! #dT&                  ,dCi_CM_dt                                     &   ! CIN            Concentration  TENDENCY, CMiPhy         [1/s]  (O)
1024! #dT&                  ,dpktCP_dt                                     &   ! Reduced Potential Temperature TENDENCY, CVAmnh        [KX/s]  (O)
1025! #dT&                  ,dqv_CP_dt                                     &   ! Specific          Humidity    TENDENCY, CVAmnh     [kg/kg/s]  (O)
1026! #dT&                  ,dqw_CP_dt                                     &   ! Cloud Droplets Concentration  TENDENCY, CVAmnh     [kg/kg/s]  (O)
1027! #dT&                  ,dqi_CP_dt                                     &   ! Cloud Crystals Concentration  TENDENCY, CVAmnh     [kg/kg/s]  (O)
1028! #dT&                  ,dpktRT_dt                                     &   ! Reduced Potential Temperature TENDENCY, radCEP        [KX/s]  (O)
1029&                  ,sst__HOST                                     &   ! Ocean     FORCING (SST)                                  [K] I
1030! #IP&                  ,sif__HOST                                     &   ! Ocean     FORCING (Sea-Ice Fraction )                    [-] I
1031! #AO&                  ,s_T__HOST                                     &   ! Ocean    COUPLING (Surface Temperat.)  n=1: Open Ocean   [-] I,NEMO
1032! #AO&                  ,Alb__HOST                                     &   ! Ocean    COUPLING (Surface Albedo   )  n=2: Sea  Ice     [-] I,NEMO
1033! #AO&                  ,dSdT2HOST                                     &   ! Ocean    COUPLING ( d(SH Flux) / dT )               [W/m2/K]   O
1034! #AO&                  ,dLdT2HOST                                     &   ! Ocean    COUPLING ( d(LH Flux) / dT )               [W/m2/K]   O
1035!dead&                  ,it0EXP,it0RUN                                 &   ! Iteration
1036&                  ,Year_H,Mon__H,Day__H,Hour_H,minu_H,sec__H     &   ! Time
1037&                  ,ixq1  ,i0x0  ,mxqq                            &   ! Domain  Dimension: x
1038&                  ,jyq1  ,j0y0  ,myqq                            &   ! Domain  Dimension: y
1039&                  ,klev   ,klevp1                                   &   ! Domain  Dimension: z
1040&                  ,mwq                                           &   ! Domain  Dimension: mosaic
1041&                  ,klon                                         &   ! Domain  Dimension: x * y                                             NEW
1042&                  ,kcolw                                         &   ! Domain  Dimension: x * y * mosaic                                    NEW
1043&                  ,m_azim                                        &   ! Mountain Mask, nb of directions taken into account       [-]         NEW
1044&                  ,IOi0SV,IOj0SV,n0pt)                               ! Indices of OUTPUT Grid Point
1045!                *******
1046
1047!======================================================================
1048! On renvoie les variables nécessaires à la dynamique de LMDZ:
1049!======================================================================
1050
1051PRINT*, 'On est sorti de PHY_MAR, on peut actualiser les variables pour LMDZ'
1052
1053d_ps(:)=0.0           ! d_ps----output-R-tendance physique de la pression au sol ! égal à zéro sur la terre (car on néglige la perte de poid de la colonne lorsqu'il pleut).
1054
1055DO k = 1, klev
1056  i=klev+1-k
1057
1058  d_u(:,k)=dua____dt(:,i) ! d_u-----output-R-tendance physique de "u" (m/s/s)
1059  d_v(:,k)=dva____dt(:,i) ! d_v-----output-R-tendance physique de "v" (m/s/s)
1060  d_t(:,k)=dpkt___dt(:,i)*(pplay(:,k)/1000) ** kap !d_t(:,:)= ! d_t-----output-R-tendance physique de "t" (K/s). Note: la pression ne change pas au cours de la physique.
1061  ! d_qx----output-R-tendance physique de "qx" (kg/kg/s):
1062  d_qx(:,k,ivap)=dqv____dt(:,i)        !  Specific           Humidity    TENDENCY, ALL Contr. [kg/kg/s]
1063  d_qx(:,k,iliq)=dqw____dt(:,i)        !  Cloud Droplets Concentration   TENDENCY, ALL Contr. [kg/kg/s]
1064  d_qx(:,k,iice)=dqi____dt(:,i)        !  Cloud Crystals Concentration   TENDENCY, ALL Contr. [kg/kg/s]
1065  d_qx(:,k,icin)=dCi____dt(:,i)        !  CIN            Concentration   TENDENCY, ALL Contr.  [1/kg/s]
1066  d_qx(:,k,isnow)=dqs____dt(:,i)       !  Snow Particles Concentration   TENDENCY, ALL Contr. [kg/kg/s]
1067  d_qx(:,k,irain)=dqr____dt(:,i)       !  Rain Drops     Concentration   TENDENCY, ALL Contr. [kg/kg/s]
1068  ! #cw dCw____dt                      !  CCN            Concentration   TENDENCY, ALL Contr.  [1/kg/s]
1069  cldfra(:,k)=CF___HOST(:,i)        ! cloud fraction
1070
1071END DO
1072
1073! Heat flux at surface (diagnostic)
1074      HsenSV_gpt(:) = -uts_SV_gpt(:) *1.e-3 *1005.    ! sensible heat flx in W/m2
1075      HLatSV_gpt(:) = -uqs_SV_gpt(:) *1.e-3 *2.5e6    ! latent heat flx in W/m2
1076
1077! Question à Hubert : tendences de TKE non nécessaires ? Réponse: la TKE est mise à jour dans la physique, on n'a pas besoin de transporter la dérivée de TKE dans la dynamique.
1078
1079! On désalloue les tableaux:
1080DEALLOCATE(MMaskHOST) ! Est-ce qu'il faut le faire ?
1081
1082PRINT*, 'Control des tendances aprés la physique de MAR'
1083
1084PRINT*, 'maxval(d_u)=',maxval(d_u)
1085PRINT*, 'maxval(d_v)=',maxval(d_v)
1086PRINT*, 'maxval(d_t)=',maxval(d_t)
1087PRINT*, 'maxval(paprs(:,1))=',maxval(paprs(:,1))
1088PRINT*, 'minval(d_u)=',minval(d_u)
1089PRINT*, 'minval(d_v)=',minval(d_v)
1090PRINT*, 'minval(d_t)=',minval(d_t)
1091PRINT*, 'minval(paprs(:,1))=',minval(paprs(:,1))
1092PRINT*,'Test snowCM'
1093PRINT*,'MAXVAL(snowCM)=',MAXVAL(snowCM)
1094PRINT*,'MAXVAL(rainCM)=',MAXVAL(rainCM)
1095PRINT*,'maxval(cldfra)=',maxval(cldfra)
1096PRINT*,'minval(cldfra)=',minval(cldfra)
1097
1098!======================================================================
1099! Ecriture des sorties netcdf
1100!======================================================================
1101
1102PRINT*, 'On ecrit les variables dans un fichier netcdf instantané'
1103! write some outputs:
1104if (debut) then ! On imprime les variables de sortie intemporelles seulement au début
1105  !call iophys_ini
1106  call histwrite_phy(nid_hist,.false.,"lonh",itau,lonh_HOST)
1107  call histwrite_phy(nid_hist,.false.,"latr",itau,latr_HOST)
1108  call histwrite_phy(nid_hist,.false.,"sst",itau,sst__HOST)
1109endif
1110if (modulo(itau,iwrite_phys)==0) then
1111  call histwrite_phy(nid_hist,.false.,"temperature",itau,t)
1112  call histwrite_phy(nid_hist,.false.,"pplay",itau,pplay)
1113  call histwrite_phy(nid_hist,.false.,"ps",itau,paprs(:,1))
1114  call histwrite_phy(nid_hist,.false.,"qx_vap",itau,qx(:,:,ivap))
1115  call histwrite_phy(nid_hist,.false.,"qx_liq",itau,qx(:,:,iliq))
1116  call histwrite_phy(nid_hist,.false.,"qx_ice",itau,qx(:,:,iice))
1117  call histwrite_phy(nid_hist,.false.,"u",itau,u)
1118  call histwrite_phy(nid_hist,.false.,"v",itau,v)
1119  call histwrite_phy(nid_hist,.false.,"d_t",itau,d_t)
1120  call histwrite_phy(nid_hist,.false.,"d_u",itau,d_u)
1121  call histwrite_phy(nid_hist,.false.,"d_v",itau,d_v)
1122  call histwrite_phy(nid_hist,.false.,"snow",itau,SnowCM)
1123  call histwrite_phy(nid_hist,.false.,"rain",itau,RainCM)
1124  call histwrite_phy(nid_hist,.false.,"snowCP",itau,snowCP)
1125  call histwrite_phy(nid_hist,.false.,"rainCP",itau,rainCP)
1126  call histwrite_phy(nid_hist,.false.,"nebulosity",itau,cldfra)
1127endif
1128
1129!XIOS
1130#ifdef CPP_XIOS
1131!$OMP MASTER
1132    !On incrémente le pas de temps XIOS
1133    !CALL wxios_update_calendar(itau)
1134
1135    !Et on écrit, avec la routine histwrite dédiée:
1136    !CALL histwrite_phy("temperature",t)
1137    !CALL histwrite_phy("u",u)
1138    !CALL histwrite_phy("v",v)
1139    !CALL histwrite_phy("ps",paprs(:,1))
1140!$OMP END MASTER
1141#endif
1142!
1143! On synchronise la fin de l'écriture du fichier à chaque pas de temps pour que le fichier soit OK même si le code plante...
1144CALL histsync(nid_hist)
1145
1146PRINT*, 'Fin de physiq.f90'
1147
1148end subroutine physiq
1149
1150
1151END MODULE physiq_mod
Note: See TracBrowser for help on using the repository browser.