4948cbac8c6b592a5ab786d734cf2ae1fff284c7
[MNH-git_open_source-lfs.git] / src / MNH / prep_ideal_case.f90
1 !MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
3 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
4 !MNH_LIC for details. version 1.
5 !-----------------------------------------------------------------
6 !--------------- special set of characters for RCS information
7 !-----------------------------------------------------------------
8 ! $Source$ $Revision$
9 !-----------------------------------------------------------------
10 !     #######################
11       PROGRAM PREP_IDEAL_CASE
12 !     #######################
13 !
14 !!****  *PREP_IDEAL_CASE* - program to write an initial FM-file 
15 !!
16 !!    PURPOSE
17 !!    -------
18 !       The purpose of this program is to prepare an initial meso-NH file
19 !     (LFIFM and DESFM files) filled with some idealized fields.    
20 !
21 !      ---- The present version can provide two types of fields:
22 !
23 !      1) CIDEAL = 'CSTN' : 3D fields derived  from a vertical profile with
24 !         ---------------   n levels of constant moist Brunt Vaisala frequency
25 !             The vertical profile is read in EXPRE file.                 
26 !             These fields can be used for model runs 
27 !
28 !      2) CIDEAL = 'RSOU' : 3D fields derived from a radiosounding.
29 !          --------------- 
30 !             The radiosounding is read in EXPRE file. 
31 !             The following kind of data  is permitted :
32 !                  YKIND = 'STANDARD'  :   Zsol, Psol, Tsol, TDsol
33 !                                         (Pressure, dd, ff) , 
34 !                                         (Pressure, T, Td)
35 !                  YKIND = 'PUVTHVMR'  : zsol, Psol, Thvsol, Rsol
36 !                                        (Pressure, U, V) , 
37 !                                        (Pressure, THv, R)
38 !                  YKIND = 'PUVTHVHU'  :  zsol, Psol, Thvsol, Husol
39 !                                         (Pressure, U, V) , 
40 !                                         (Pressure, THv, Hu)
41 !                  YKIND = 'ZUVTHVHU'  :  zsol, Psol, Thvsol, Husol
42 !                                         (height, U, V) , 
43 !                                         (height, THv, Hu)
44 !                  YKIND = 'ZUVTHVMR'  :  zsol, Psol, Thvsol, Rsol
45 !                                         (height, U, V) , 
46 !                                         (height, THv, R)
47 !                  YKIND = 'PUVTHDMR'  : zsol, Psol, Thdsol, Rsol
48 !                                         (Pressure, U, V) , 
49 !                                         (Pressure, THd, R)
50 !                  YKIND = 'PUVTHDHU'  : zsol, Psol, Thdsol, Husol
51 !                                         (Pressure, U, V) , 
52 !                                         (Pressure, THd, Hu)
53 !                  YKIND = 'ZUVTHDMR'  :  zsol, Psol, Thdsol, Rsol
54 !                                         (height, U, V) , 
55 !                                         (height, THd, R)
56 !                  YKIND = 'ZUVTHLMR'  :  zsol, Psol, Thdsol, Rsol
57 !                                         (height, U, V) , 
58 !                                         (height, THl, Rt)
59 !
60 !             These fields can be used for model runs 
61 !
62 !      Cases (1) and (2) can be balanced
63 !      (geostrophic, hydrostatic  and anelastic balances) if desired.
64 !
65 !      ---- The orography can be flat (YZS='FLAT'), but also 
66 !      sine-shaped (YZS='SINE') or  bell-shaped (YZS='BELL')
67 !
68 !      ---- The U(z)  profile given in the RSOU and CSTN cases can
69 !      be multiplied (CUFUN="Y*Z") by a function of y (function FUNUY)  
70 !      The V(z) profile  given in the RSOU and CSTN cases can
71 !      be multiplied (CVFUN="X*Z") by a function of x (function FUNVX). 
72 !      If it is not the case, i.e. U(y,z)=U(z) then CUFUN="ZZZ" and 
73 !      CVFUN="ZZZ" for V(y,z)=V(z). Instead of these separable forms,
74 !      non-separables functions FUNUYZ (CUFUN="Y,Z")  and FUNVXZ (CVFUN="X,Z") 
75 !      can be used to specify the wind components.
76 !
77 !!**  METHOD
78 !!    ------
79 !!      The directives and data to perform the preparation of the initial FM
80 !!    file are stored in EXPRE file. This file is composed  of two parts : 
81 !!          - a namelists-format  part which is present in all cases
82 !!          - a free-format  part which contains data in cases 
83 !!       of discretised orography (CZS='DATA')
84 !!       of radiosounding (CIDEAL='RSOU') or Nv=cste  profile (CIDEAL='CSTN')
85 !!       of forced version (LFORCING=.TRUE.)
86 !!    
87 !!
88 !!      The following  PREP_IDEAL_CASE program  :
89 !!
90 !!             - initializes physical constants by calling INI_CST 
91 !!
92 !!             - sets default values for global variables which will be 
93 !!     written  in DESFM file and for variables in EXPRE file (namelists part)
94 !!     which will be written in LFIFM file.    
95 !!
96 !!             - reads the namelists part of EXPRE file which gives 
97 !!     informations about the preinitialization to perform,
98 !!
99 !!             - allocates memory for arrays, 
100 !!
101 !!             - initializes fields depending on the 
102 !!              directives  (CIDEAL in namelist NAM_CONF_PRE) :
103 !!  
104 !!                * grid variables : 
105 !!                  The gridpoints are regularly spaced by XDELTAX, XDELTAY.
106 !!               The grid is stretched along the z direction, the mesh varies 
107 !!               from XDZGRD near the ground to XDZTOP near the top and the 
108 !!               weigthing function is a TANH function characterized by its 
109 !!               center and width above and under this center
110 !!                  The orography is initialized following the kind of orography
111 !!               (YZS in namelist NAM_CONF_PRE) and the degrees of freedom :
112 !!                     sine-shape ---> ZHMAX, IEXPX,IEXPY
113 !!                     bell-shape ---> ZHMAX, ZAX,ZAY,IIZS,IJZS
114 !!                  The horizontal grid variables are initialized following
115 !!                the kind of geometry (LCARTESIAN in namelist NAM_CONF_PRE) 
116 !!                and the grid parameters XLAT0,XLON0,XBETA in both geometries
117 !!                and XRPK,XLONORI,XLATORI  in conformal projection.
118 !!                  In the  case of initialization from a radiosounding, the
119 !!                date and time is read in free-part of the EXPRE file. In other
120 !!                cases year, month and day are set to NUNDEF and time to 0.
121 !!
122 !!               * prognostic fields : 
123 !!
124 !!                     U,V,W, Theta and r. are first determined. They are
125 !!                multiplied by rhoj after the anelastic reference state 
126 !!                computation.
127 !!                     For the CSTN and RSOU cases, the determination of 
128 !!                Theta and rv is performed  respectively by SET_RSOU
129 !!                and by SET_CSTN which call the common routine SET_MASS. 
130 !!                These three routines have  the following actions :
131 !!          ---   The input vertical profile   is converted in 
132 !!                variables (U,V,thetav,r) and  interpolated
133 !!                on a mixed grid (with VERT_COORD) as in PREP_REAL_CASE 
134 !!          ---   A variation of the u-wind component( x-model axis component) 
135 !!                 is possible in y direction, a variation of the v-wind component 
136 !!                (y-model axis component) is possible in x direction.
137 !!          ---   Thetav could be computed with thermal wind balance
138 !!                (LGEOSBAL=.TRUE. with call of SET_GEOSBAL)                 
139 !!          ---   The mass fields (theta and r ) and the wind components are 
140 !!                then interpolated on the model grid with orography  as in
141 !!                PREP_REAL_CASE with the option LSHIFT                
142 !!          ---   An  anelastic correction is  applied in PRESSURE_IN_PREP in
143 !!                the case of non-vanishing orography.    
144 !!            
145 !!               * anelastic reference state variables :
146 !!
147 !!                   1D reference state : 
148 !!                     RSOU and CSTN cases : rhorefz and thvrefz are computed 
149 !!                         by   SET_REFZ (called by SET_MASS).
150 !!                         They are deduced from thetav and r on the model grid
151 !!                         without orography.
152 !!                   The 3D reference state is  computed by SET_REF   
153 !!            
154 !!               * The total mass of dry air is computed by TOTAL_DMASS              
155 !!
156 !!             - writes the DESFM file, 
157 !!
158 !!             - writes the LFIFM file . 
159 !!
160 !!    EXTERNAL
161 !!    --------
162 !!      DEFAULT_DESFM : to set default values for variables which can be 
163 !!                      contained in DESFM file
164 !!      DEFAULT_EXPRE : to  set default values for other global variables 
165 !!                      which can be contained in namelist-part of EXPRE file
166 !!      Module MODE_GRIDPROJ : contains conformal projection routines
167 !!           SM_GRIDPROJ   : to compute some grid variables, in
168 !!                           case of conformal projection.
169 !!      Module MODE_GRIDCART : contains cartesian geometry routines
170 !!           SM_GRIDCART   : to compute some grid variables, in
171 !!                           case of cartesian geometry.
172 !!      SET_RSOU      : to initialize mass fields from a radiosounding
173 !!      SET_CSTN      : to initialize mass fields from a vertical profile of 
174 !!                      n layers of Nv=cste 
175 !!      SET_REF       : to compute  rhoJ 
176 !!      RESSURE_IN_PREP : to apply an anelastic correction in the case of
177 !!                        non-vanishing orography 
178 !!      FMOPEN        : to open a FM-file (DESFM + LFIFM)
179 !!      WRITE_DESFM   : to write the  DESFM file
180 !!      WRI_LFIFM     : to write the   LFIFM file  
181 !!      FMCLOS        : to close a FM-file (DESFM + LFIFM)
182 !!
183 !!      MXM,MYM,MZM   : Shuman operators
184 !!      WGUESS        : to compute W with the continuity equation from 
185 !!                      the U,V values 
186 !!
187 !!
188 !!      
189 !!    IMPLICIT ARGUMENTS
190 !!    ------------------
191 !!      Module MODD_PARAMETERS : contains parameters
192 !!      Module MODD_DIM1      : contains dimensions 
193 !!      Module MODD_CONF       : contains  configuration variables for 
194 !!                                all models
195 !!      Module MODD_CST        : contains physical constants
196 !!      Module MODD_GRID       : contains grid variables  for all models
197 !!      Module MODD_GRID1     : contains grid variables
198 !!      Module MODD_TIME      : contains time variables for all models  
199 !!      Module MODD_TIME1     : contains time variables  
200 !!      Module MODD_REF        : contains reference state variables for
201 !!                               all models
202 !!      Module MODD_REF1      : contains reference state variables 
203 !!      Module MODD_LUNIT      : contains variables which concern names
204 !!                            and logical unit numbers of files  for all models
205 !!      Module MODD_FIELD1    : contains prognostics  variables
206 !!      Module MODD_GR_FIELD1 : contains the surface prognostic variables 
207 !!      Module MODD_LSFIELD1    : contains Larger Scale fields
208 !!      Module MODD_DYN1        : contains dynamic control variables for model 1
209 !!      Module MODD_LBC1        : contains lbc control variables for model 1
210 !!
211 !!
212 !!      Module MODN_CONF1    : contains  configuration variables for model 1
213 !!                               and the NAMELIST list
214 !!      Module MODN_LUNIT1    : contains variables which concern names
215 !!                               and logical unit numbers of files and 
216 !!                               the NAMELIST list
217 !!
218 !!
219 !!    REFERENCE
220 !!    ---------
221 !!      Book2 of MESO-NH documentation (program PREP_IDEAL_CASE)
222 !!    
223 !!    AUTHOR
224 !!    ------
225 !!      V. Ducrocq   *Meteo France*
226 !!
227 !!    MODIFICATIONS
228 !!    -------------
229 !!      Original                            05/05/94
230 !!      updated                V. Ducrocq   27/06/94   
231 !!      updated                P.M.         27/07/94
232 !!      updated                V. Ducrocq   23/08/94 
233 !!      updated                V. Ducrocq   01/09/94 
234 !!      namelist changes       J. Stein     26/10/94 
235 !!      namelist changes       J. Stein     04/11/94 
236 !!      remove the second step of the geostrophic balance 14/11/94 (J.Stein)
237 !!      add grid stretching in the z direction + Larger scale fields +
238 !!      cleaning                                           6/12/94 (J.Stein) 
239 !!      periodize the orography and the grid sizes in the periodic case
240 !!                                                        19/12/94 (J.Stein) 
241 !!      correct a bug in the Larger Scale Fields initialization
242 !!                                                        19/12/94 (J.Stein) 
243 !!      add the vertical grid stretching                  02/01/95 (J. Stein)
244 !!      Total mass of dry air computation                 02/01/95 (J.P.Lafore) 
245 !!      add the 1D switch                                 13/01/95 (J. Stein)
246 !!      enforce a regular vertical grid if desired        18/01/95 (J. Stein)
247 !!      add the tdtcur initialization                     26/01/95 (J. Stein)
248 !!      bug in the test of the type of RS localization    25/02/95 (J. Stein)
249 !!      remove R from the historical variables            16/03/95 (J. Stein)
250 !!      error on the grid stretching                      30/06/95 (J. Stein)
251 !!      add the soil fields                               01/09/95 (S.Belair)
252 !!      change the streching function  and the wind guess
253 !!        (J. Stein and V.Masson)                         21/09/95 
254 !!      reset to FALSE LUSERC,..,LUSERH                   12/12/95 (J. Stein)
255 !!      enforce the RS localization in 1D and 2D config.
256 !!      + add the 'TSZ0' option for the soil variables    28/01/96 (J. Stein)
257 !!      initialization of domain from center point        31/01/96 (V. Masson)
258 !!      add the constant file reading                     05/02/96 (J. Stein)
259 !!      enter vertical model levels values                20/10/95 (T.Montmerle)
260 !!      add LFORCING option                               19/02/96 (K. Suhre)
261 !!      modify structure of NAM_CONF_PRE                  20/02/96 (J.-P. Pinty)
262 !!      default of the domain center when use of pgd file 12/03/96 (V. Masson)
263 !!      change the surface initialization                 20/03/96 ( Stein,
264 !!                                                    Bougeault, Kastendeutsch )
265 !!      change the DEFAULT_DESFMN CALL                    17/04/96 ( Lafore )
266 !!      set the STORAGE_TYPE to 'TT' (a single instant)   30/04/96 (Stein, 
267 !!                                                    Jabouille)
268 !!      new wguess to spread  the divergence              15/05/96 (Stein)
269 !!      set LTHINSHELL to TRUE + return to the old wguess 29/08/96 (Stein)
270 !!      MY_NAME and DAD_NAME writing for nesting          30/07/96 (Lafore)
271 !!      MY_NAME and DAD_NAME reading in pgd file          26/09/96 (Masson)
272 !!       and reading of pgd grid in a new routine
273 !!      XXHAT and XYHAT are set to 0. at origine point    02/10/96 (Masson)
274 !!      add LTHINSHELL in namelist NAM_CONF_PRE           08/10/96 (Masson)
275 !!      restores use of TS and T2                         26/11/96 (Masson)
276 !!      value  XUNDEF for soil and vegetation fields on sea 27/11/96 (Masson)
277 !!      use of HUG and HU2 in both ISBA and TSZ0 cases    04/12/96 (Masson)
278 !!      add initialization of chemical variables          06/08/96 (K. Suhre)
279 !!      add MANUAL option for the terrain elevation       12/12/96 (J.-P. Pinty)
280 !!      set DATA instead of MANUAL for the terrain
281 !!      elevation option
282 !!      add new anelastic equations' systems              29/06/97 (Stein)
283 !!      split mode_lfifm_pgd                              29/07/97 (Masson)
284 !!      add directional z0 and subgrid scale orography    31/07/97 (Masson)
285 !!      separates surface treatment in PREP_IDEAL_SURF    15/03/99 (Masson)
286 !!      new PGD fields allocations                        15/03/99 (Masson)
287 !!      iterative call to pressure solver                 15/03/99 (Masson)
288 !!      removes TSZ0 case                                 04/01/00 (Masson)
289 !!      parallelization                                   18/06/00 (Pinty)
290 !!      adaptation for patch approach                     02/07/00 (Solmon/Masson)
291 !!      bug in W LB field on Y direction                  05/03/01 (Stein)
292 !!      add module MODD_NSV for NSV variable              01/02/01 (D. Gazen) 
293 !!      allow namelists in different orders               15/10/01 (I. Mallet)
294 !!      allow LUSERC and LUSERI in 1D configuration       05/06/02 (P. Jabouille)
295 !!      add  ZUVTHLMR case (move in set_rsou latter)      05/12/02 Jabouille/Masson
296 !!      move LHORELAX_SV (after INI_NSV)                  30/04/04 (Pinty)
297 !!      Correction Parallel bug IBEG & IDEND  evalution   13/11/08 J.Escobar
298 !!      add the option LSHIFT for interpolation of        26/10/10 (G.Tanguy)
299 !!      correction for XHAT & parallelizarion of ZSDATA   23/09/11 J.Escobar
300 !!      the vertical profile (as in PREP_REAL_CASE)
301 !!      add use MODI of SURFEX routines                   10/10/111 J.Escobar
302 !!
303 !!      For 2D modeling: 
304 !!      Initialization of ADVFRC profiles (SET_ADVFRC)    06/2010 (P.Peyrille)
305 !!      when LDUMMY(2)=T in PRE_IDEA1.nam 
306 !!      USE MODDB_ADVFRC_n for grid-nesting               02*2012 (M. Tomasini)
307 !!      LBOUSS in MODD_REF                                07/2013 (C.Lac)
308 !!      Correction for ZS in PGD file                     04/2014 (G. TANGUY)
309 !!      Bug : remove NC WRITE_HGRID                       05/2014 (S. Bielli via J.Escobar )
310 !!      BUG if ZFRC and ZFRC_ADV or ZFRC_REL are used together  11/2014 (G. Delautier)
311 !!      Bug : detected with cray compiler ,
312 !!                  missing '&' in continuation string  3/12/2014 J.Escobar
313 !!      J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
314 !-------------------------------------------------------------------------------
315 !
316 !*       0.   DECLARATIONS
317 !             ------------
318 !
319 USE MODD_PARAMETERS       ! Declarative modules
320 USE MODD_DIM_n
321 USE MODD_CONF
322 USE MODD_CST
323 USE MODD_GRID         
324 USE MODD_GRID_n
325 USE MODD_METRICS_n
326 USE MODD_PGDDIM
327 USE MODD_PGDGRID
328 USE MODD_TIME
329 USE MODD_TIME_n
330 USE MODD_REF
331 USE MODD_REF_n
332 USE MODD_LUNIT
333 USE MODD_FIELD_n
334 USE MODD_DYN_n
335 USE MODD_LBC_n
336 USE MODD_LSFIELD_n
337 USE MODD_PARAM_n
338 USE MODD_CH_MNHC_n, ONLY:  LUSECHEM, LUSECHAQ, LUSECHIC, LCH_PH,  &
339                            LCH_INIT_FIELD, CCHEM_INPUT_FILE 
340 USE MODD_CH_AEROSOL,ONLY:  LORILAM, CORGANIC, LVARSIGI, LVARSIGJ, LINITPM, XINIRADIUSI, &
341                            XINIRADIUSJ, XINISIGI, XINISIGJ, XN0IMIN, XN0JMIN, CRGUNIT
342 USE MODD_DUST,      ONLY:  LDUST, NMODE_DST, CRGUNITD, XINISIG, XINIRADIUS, XN0MIN 
343 USE MODD_SALT,      ONLY:  LSALT, NMODE_SLT, CRGUNITS, XINISIG_SLT, XINIRADIUS_SLT, XN0MIN_SLT
344 USE MODD_VAR_ll,    ONLY:  NPROC
345 USE MODD_LUNIT_n
346 USE MODD_CONF_n
347 USE MODD_NSV,      ONLY : NSV,NSV_CHEM,           &
348                           NSV_DSTEND, NSV_DSTBEG
349 !
350 USE MODN_BLANK
351 !
352 USE MODE_THERMO
353 USE MODE_POS
354 USE MODE_GRIDCART         ! Executive modules
355 USE MODE_GRIDPROJ
356 USE MODE_FM
357 USE MODE_FMREAD
358 USE MODE_IO_ll
359 USE MODE_ll
360 USE MODD_ARGSLIST_ll, ONLY : LIST_ll
361 USE MODE_MODELN_HANDLER
362 !
363 USE MODI_DEFAULT_DESFM_n    ! Interface modules
364 USE MODI_DEFAULT_EXPRE
365 USE MODI_READ_HGRID
366 USE MODI_SHUMAN
367 USE MODI_SET_RSOU
368 USE MODI_SET_CSTN
369 USE MODI_SET_FRC
370 USE MODI_PRESSURE_IN_PREP
371 USE MODI_WRITE_DESFM_n
372 USE MODI_WRITE_LFIFM_n
373 USE MODI_METRICS
374 USE MODI_UPDATE_METRICS
375 USE MODI_SET_REF
376 USE MODI_SET_PERTURB
377 USE MODI_TOTAL_DMASS
378 USE MODI_WGUESS
379 USE MODI_CH_INIT_SCHEME_n
380 USE MODI_CH_INIT_FIELD_n
381 USE MODI_GATHER_ll
382 USE MODI_INI_NSV
383 USE MODI_READ_PRE_IDEA_NAM_n
384 USE MODI_CH_AER_INIT_SOA
385 USE MODI_ZSMT_PIC
386 USE MODI_ZSMT_PGD
387 USE MODI_READ_VER_GRID
388 USE MODI_READ_ALL_NAMELISTS
389 USE MODI_GOTO_SURFEX
390 USE MODI_PGD_GRID_SURF_ATM
391 USE MODI_SPLIT_GRID
392 USE MODI_PGD_SURF_ATM
393 USE MODI_ICE_ADJUST_BIS
394 USE MODI_WRITE_PGD_SURF_ATM_n
395 USE MODI_PREP_SURF_MNH
396 USE MODI_ALLOC_SURFEX
397 !
398 !JUAN
399 USE MODE_SPLITTINGZ_ll
400 USE MODD_SUB_MODEL_n
401 USE MODE_MNH_TIMING
402 USE MODN_CONFZ
403 !JUAN
404 #ifdef MNH_NCWRIT
405 USE MODN_NCOUT
406 USE MODE_UTIL
407 #endif
408 USE MODN_CONFIO
409 USE MODI_TH_R_FROM_THL_RT_3D
410 !
411 USE MODI_VERSION
412 USE MODI_INIT_PGD_SURF_ATM
413 USE MODI_WRITE_SURF_ATM_N
414 USE MODI_DEALLOC_SURF_ATM_N
415 USE MODI_DEALLOC_SURFEX
416 ! Modif ADVFRC
417 USE MODD_2D_FRC
418 USE MODD_ADVFRC_n     ! Modif for grid-nesting
419 USE MODI_SETADVFRC
420 USE MODD_RELFRC_n     ! Modif for grid-nesting
421 USE MODI_SET_RELFRC
422 !
423 USE MODI_INI_CST
424 USE MODI_INI_NEB
425 USE MODE_FMWRIT
426 USE MODI_WRITE_HGRID
427 USE MODD_MPIF
428 USE MODD_VAR_ll
429 !
430 USE MODE_MPPDB
431 !
432 IMPLICIT NONE
433 !
434 !*       0.1  Declarations of global variables not declared in the modules
435 !
436 REAL, DIMENSION(:,:,:), ALLOCATABLE :: XJ ! Jacobian
437 REAL :: XLATCEN=XUNDEF, XLONCEN=XUNDEF ! latitude and longitude of the center of
438                                      ! the domain for initialization. This 
439                                      ! point is vertical vorticity point
440                                      !          ------------------------
441 REAL :: XDELTAX=0.5E4, XDELTAY=0.5E4 ! horizontal mesh lengths  
442                                      !  used to determine  XXHAT,XYHAT
443 !
444 INTEGER :: NLUPRE,NLUOUT           ! Logical unit numbers for EXPRE file
445                                    ! and for output_listing file
446 INTEGER :: NRESP                   ! return code in FM routines
447 INTEGER :: NTYPE                   ! type of file (cpio or not)
448 INTEGER :: NNPRAR                  ! number of articles predicted  in
449                                    !  the LFIFM file
450 INTEGER :: NNINAR                  ! number of articles  present in
451                                    !  the LFIFM file
452 LOGICAL :: GFOUND                  ! Return code when searching namelist
453 !
454 INTEGER :: JLOOP,JILOOP,JJLOOP     ! Loop indexes
455 !
456 INTEGER :: NIB,NJB,NKB             ! Begining useful area  in x,y,z directions
457 INTEGER :: NIE,NJE                 ! Ending useful area  in x,y directions
458 INTEGER :: NIU,NJU,NKU             ! Upper bounds in x,y,z directions
459 CHARACTER (LEN=32) ::  CEXPRE            ! name of the EXPRE file
460 CHARACTER (LEN=32) :: CDESFM             ! Name of DESFM file 
461 CHARACTER(LEN=4)   :: CIDEAL ='CSTN'     ! kind of idealized fields
462                                          ! 'CSTN' : Nv=cste case 
463                                          ! 'RSOU' : radiosounding case
464 CHARACTER(LEN=4)   :: CZS    ='FLAT'     ! orography selector
465                                          ! 'FLAT' : zero orography
466                                          ! 'SINE' : sine-shaped orography 
467                                          ! 'BELL' : bell-shaped orography 
468 REAL    :: XHMAX=XUNDEF            ! Maximum height for orography
469 REAL    :: NEXPX=3,NEXPY=1         ! Exponents for  orography in case of CZS='SINE'
470 REAL    :: XAX= 1.E4, XAY=1.E4     ! Widths for orography in case CZS='BELL'
471                                    ! along x and y 
472 INTEGER :: NIZS = 5, NJZS = 5      ! Localization of the center in 
473                                    ! case CZS ='BELL' 
474 !
475 !*       0.1.1 Declarations of local variables for N=cste and 
476 !              radiosounding cases :
477 !
478 INTEGER            :: NYEAR,NMONTH,NDAY ! year, month and day in EXPRE file
479 REAL               :: XTIME             ! time in EXPRE file
480 LOGICAL            :: LPERTURB =.FALSE. ! Logical to add a perturbation to 
481                                         ! a basic state 
482 LOGICAL            :: LGEOSBAL =.FALSE. ! Logical to satisfy the geostrophic
483                                         ! balance
484                                         ! .TRUE. for geostrophic balance
485                                         ! .FALSE. to ignore this balance
486 LOGICAL            :: LPV_PERT =.FALSE. ! Logical to add a PV pertubation
487 LOGICAL            :: LRMV_BL  =.FALSE. ! Logical to remove the boundary layer
488                                         ! before PV inversion
489 LOGICAL            :: LSHIFT   =.FALSE.  ! flag to perform vertical shift or not.        
490 CHARACTER(LEN=3)   :: CFUNU ='ZZZ'      ! CHARACTER STRING for variation of
491                                         ! U in y direction
492                                         ! 'ZZZ'  : U = U(Z)
493                                         ! 'Y*Z'  : U = F(Y) * U(Z)
494                                         ! 'Y,Z'  : U = G(Y,Z)
495 CHARACTER(LEN=3)   :: CFUNV ='ZZZ'      ! CHARACTER STRING for variation of
496                                         ! V in x direction
497                                         ! 'ZZZ'  : V = V(Z)
498                                         ! 'Y*Z'  : V = F(X) * V(Z)
499                                         ! 'Y,Z'  : V = G(X,Z)
500 CHARACTER(LEN=6)   :: CTYPELOC='IJGRID' ! Type of informations  used to give the
501                                         ! localization of vertical profile
502                                         ! 'IJGRID'  for (i,j) point  on index space
503                                         ! 'XYHATM' for (x,y) coordinates on
504                                         !  conformal or cartesian plane
505                                         ! 'LATLON' for (latitude,longitude) on
506                                         !   spherical earth  
507 REAL               :: XLATLOC= 45., XLONLOC=0.
508                                         ! Latitude and longitude of the vertical
509                                         ! profile localization  (used in case 
510                                         ! CTYPELOC='LATLON') 
511 REAL               :: XXHATLOC=2.E4, XYHATLOC=2.E4 
512                                         ! (x,y) of the vertical profile
513                                         ! localization  (used in cases 
514                                         ! CTYPELOC='LATLON' and 'XYHATM') 
515 INTEGER, DIMENSION(1) :: NILOC=4, NJLOC=4 
516                                         ! (i,j) of the vertical profile
517                                         ! localization 
518 !
519 !
520 REAL,DIMENSION(:,:,:),ALLOCATABLE   :: XCORIOZ ! Coriolis parameter (this
521                                                  ! is exceptionnaly a 3D array
522                                                  ! for computing needs)
523 !
524 !
525 !*       0.1.2 Declarations of local variables used when a PhysioGraphic Data
526 !              file is used :
527 !
528 INTEGER             :: JSV                      ! loop index on scalar var.
529 CHARACTER(LEN=28)   :: CPGD_FILE=' '            ! Physio-Graphic Data file name
530 LOGICAL  :: LREAD_ZS = .TRUE.,                & ! switch to use orography 
531                                                 ! coming from the PGD file
532             LREAD_GROUND_PARAM = .TRUE.         ! switch to use soil parameters
533                                                 ! useful for the soil scheme
534                                                 ! coming from the PGD file
535
536 INTEGER           :: NSLEVE   =12         ! number of iteration for smooth orography
537 REAL              :: XSMOOTH_ZS = XUNDEF  ! optional uniform smooth orography for SLEVE coordinate
538 CHARACTER(LEN=28) :: YPGD_NAME, YPGD_DAD_NAME   ! general information
539 CHARACTER(LEN=8)  :: YKIND                      ! Kind of radiosounding data
540 CHARACTER(LEN=2)  :: YPGD_TYPE
541 !
542 INTEGER           :: IINFO_ll                   ! return code of // routines
543 TYPE(LIST_ll), POINTER :: TZ_FIELDS_ll           ! list of metric coefficient fields
544 !
545 INTEGER :: IISIZEXF,IJSIZEXF,IISIZEXFU,IJSIZEXFU     ! dimensions of the
546 INTEGER :: IISIZEX4,IJSIZEX4,IISIZEX2,IJSIZEX2       ! West-east LB arrays
547 INTEGER :: IISIZEYF,IJSIZEYF,IISIZEYFV,IJSIZEYFV     ! dimensions of the
548 INTEGER :: IISIZEY4,IJSIZEY4,IISIZEY2,IJSIZEY2       ! North-south LB arrays
549 INTEGER :: IBEG,IEND,IXOR,IXDIM,IYOR,IYDIM,ILBX,ILBY
550 REAL, DIMENSION(:),   ALLOCATABLE :: ZXHAT_ll, ZYHAT_ll
551 !
552 REAL, DIMENSION(:,:,:), ALLOCATABLE ::ZTHL,ZT,ZRT,ZFRAC_ICE,&
553                                       ZEXN,ZLVOCPEXN,ZLSOCPEXN,ZCPH, &
554                                       ZRSATW, ZRSATI
555                                  ! variables for adjustement
556 INTEGER             :: ILENCH, IGRID, IRESP
557 CHARACTER (LEN=100) :: YCOMMENT
558 REAL                :: ZDIST
559 !
560 !JUAN TIMING
561 REAL*8,DIMENSION(2)         :: ZTIME1,ZTIME2,ZEND,ZTOT
562 CHARACTER                 :: YMI
563 INTEGER                   :: IMI
564 INTEGER::JK                                 
565 !JUAN TIMING
566 !
567 REAL, DIMENSION(:),   ALLOCATABLE :: ZZS_ll
568 INTEGER                           :: IJ 
569 INTEGER           :: NZSFILTER=1          ! number of iteration for filter for fine   orography
570 LOGICAL           :: LHSLOP=.FALSE.       ! filtering of slopes higher than XHSLOP   
571 REAL              :: XHSLOP=1.2           ! if LHSLOP filtering of slopes higher than XHSLOP
572
573 !
574 REAL              :: ZZS_MAX, ZZS_MAX_ll
575 INTEGER           :: IJPHEXT
576 !
577 !
578 !*       0.2  Namelist declarations
579 !
580 NAMELIST/NAM_CONF_PRE/ LTHINSHELL,LCARTESIAN,    &! Declarations in MODD_CONF
581                        LPACK,                    &!
582                        NVERB,CIDEAL,CZS,         &!+global variables initialized
583                        LBOUSS,LPERTURB,LPV_PERT, &! at their declarations
584                        LRMV_BL,LFORCING,CEQNSYS, &! at their declarations
585                        LSHIFT,L2D_ADV_FRC,L2D_REL_FRC, &
586                        NHALO , JPHEXT
587 NAMELIST/NAM_GRID_PRE/ XLON0,XLAT0,            & ! Declarations in MODD_GRID
588                        XBETA,XRPK,             & 
589                        XLONORI,XLATORI
590 NAMELIST/NAM_GRIDH_PRE/ XLATCEN,XLONCEN,       & ! local variables  initialized
591                  XDELTAX,XDELTAY,              & ! at their declarations
592                  XHMAX,NEXPX,NEXPY,            &
593                  XAX,XAY,NIZS,NJZS
594 NAMELIST/NAM_VPROF_PRE/LGEOSBAL, CFUNU,CFUNV,   &! global variables initialized
595                      CTYPELOC,XLATLOC,XLONLOC,  &!  at their declarations
596                      XXHATLOC,XYHATLOC,NILOC,NJLOC
597 NAMELIST/NAM_REAL_PGD/CPGD_FILE,                 & ! Physio-Graphic Data file
598                                                    !  name
599                       LREAD_ZS,                  & ! switch to use orography 
600                                                    ! coming from the PGD file
601                       LREAD_GROUND_PARAM
602 NAMELIST/NAM_SLEVE/NSLEVE, XSMOOTH_ZS
603 !
604 !*       0.3  Auxillary Namelist declarations
605 !
606 NAMELIST/NAM_AERO_PRE/ LORILAM, LINITPM, XINIRADIUSI, XINIRADIUSJ, &
607                        XINISIGI, XINISIGJ, XN0IMIN, XN0JMIN, CRGUNIT, &
608                        LDUST, LSALT, CRGUNITD, CRGUNITS,&
609                        NMODE_DST, XINISIG, XINIRADIUS, XN0MIN,&
610                        XINISIG_SLT, XINIRADIUS_SLT, XN0MIN_SLT, &
611                        NMODE_SLT
612 !
613 !-------------------------------------------------------------------------------
614 !
615 !*       0.    PROLOGUE
616 !              --------
617 CALL MPPDB_INIT()
618 !
619 CALL GOTO_MODEL(1)
620 !
621 CALL INITIO_ll()
622 NULLIFY(TZ_FIELDS_ll)
623 CALL VERSION
624 CPROGRAM='IDEAL '
625 !
626 !JUAN TIMING
627   XT_START     = 0.0
628   XT_STORE     = 0.0
629 !
630   CALL SECOND_MNH2(ZEND)
631 !
632 !JUAN TIMING
633 !
634 !*       1.    INITIALIZE PHYSICAL CONSTANTS :         
635 !              ------------------------------
636 !
637 NVERB = 5
638 CALL INI_CST
639 CALL INI_NEB
640 !
641 !-------------------------------------------------------------------------------
642 !
643 !
644 !*       2.    SET DEFAULT VALUES  :  
645 !              --------------------
646 !
647 !
648 !*       2.1  For variables in DESFM file
649 !
650 CALL DEFAULT_DESFM_n(1)
651 !
652 CSURF = "NONE"
653 !
654 !
655 !*       2.2  For other global variables in EXPRE file
656 !
657 CALL DEFAULT_EXPRE
658 !-------------------------------------------------------------------------------
659 !
660 !*       3.    READ THE EXPRE FILE :  
661 !              --------------------
662 !
663 !*       3.1   initialize logical unit numbers (EXPRE and output-listing files)
664 !              and open these files :
665
666
667 CLUOUT  = 'OUTPUT_LISTING1'
668 CLUOUT0 = CLUOUT
669 CEXPRE  = 'PRE_IDEA1.nam'
670 CALL OPEN_ll(UNIT=NLUOUT,FILE=CLUOUT,IOSTAT=NRESP,FORM='FORMATTED',ACTION='WRITE', &
671      MODE=GLOBAL)
672 CALL OPEN_ll(UNIT=NLUPRE,FILE=CEXPRE,IOSTAT=NRESP,ACTION='READ', &
673      DELIM='QUOTE',MODE=GLOBAL)    
674 !
675 !*       3.2   read in NLUPRE the namelist informations
676 !
677 WRITE(NLUOUT,FMT=*) 'attempt to read ',TRIM(CEXPRE),' file'
678 CALL POSNAM(NLUPRE,'NAM_REAL_PGD',GFOUND,NLUOUT)
679 IF (GFOUND) READ(UNIT=NLUPRE,NML=NAM_REAL_PGD)
680 !
681 !
682 CALL POSNAM(NLUPRE,'NAM_CONF_PRE',GFOUND,NLUOUT)
683 IF (GFOUND) READ(UNIT=NLUPRE,NML=NAM_CONF_PRE)
684 !JUANZ
685 CALL POSNAM(NLUPRE,'NAM_CONFZ',GFOUND,NLUOUT)
686 IF (GFOUND) READ(UNIT=NLUPRE,NML=NAM_CONFZ)
687 !JUANZ
688 #ifdef MNH_NCWRIT
689 CALL POSNAM(NLUPRE,'NAM_NCOUT',GFOUND,NLUOUT)
690 IF (GFOUND) READ(UNIT=NLUPRE,NML=NAM_NCOUT)
691 #endif
692 CALL POSNAM(NLUPRE,'NAM_CONFIO',GFOUND,NLUOUT)
693 IF (GFOUND) READ(UNIT=NLUPRE,NML=NAM_CONFIO)
694 CALL SET_CONFIO_ll(LCDF4, LLFIOUT, LLFIREAD)
695 CALL POSNAM(NLUPRE,'NAM_GRID_PRE',GFOUND,NLUOUT)
696 IF (GFOUND) READ(UNIT=NLUPRE,NML=NAM_GRID_PRE)
697 CALL POSNAM(NLUPRE,'NAM_GRIDH_PRE',GFOUND,NLUOUT)
698 IF (GFOUND) READ(UNIT=NLUPRE,NML=NAM_GRIDH_PRE)
699 CALL POSNAM(NLUPRE,'NAM_VPROF_PRE',GFOUND,NLUOUT)
700 IF (GFOUND) READ(UNIT=NLUPRE,NML=NAM_VPROF_PRE)
701 CALL POSNAM(NLUPRE,'NAM_BLANK',GFOUND,NLUOUT)
702 IF (GFOUND) READ(UNIT=NLUPRE,NML=NAM_BLANK)
703 CALL READ_PRE_IDEA_NAM_n(NLUPRE,NLUOUT)
704 CALL POSNAM(NLUPRE,'NAM_AERO_PRE',GFOUND,NLUOUT)
705 IF (GFOUND) READ(UNIT=NLUPRE,NML=NAM_AERO_PRE)
706 !
707 IF( LEN_TRIM(CPGD_FILE) /= 0 ) THEN 
708   ! open the PGD_FILE
709   CALL FMOPEN_ll(CPGD_FILE,'READ',CLUOUT,NNPRAR,2,NVERB,NNINAR,NRESP)
710   ! read the grid in the PGD file
711   CALL FMREAD(CPGD_FILE,'IMAX',CLUOUT,'--',NIMAX,IGRID,ILENCH,YCOMMENT,IRESP)
712   CALL FMREAD(CPGD_FILE,'JMAX',CLUOUT,'--',NJMAX,IGRID,ILENCH,YCOMMENT,IRESP)
713   CALL FMREAD(CPGD_FILE,'JPHEXT',CLUOUT,'--',IJPHEXT,IGRID,ILENCH,YCOMMENT,IRESP)
714
715   IF ( CPGD_FILE /= CINIFILEPGD) THEN
716      WRITE(NLUOUT,FMT=*) ' WARNING : in PRE_IDEA1.nam, in NAM_LUNITn you&
717           & have CINIFILEPGD= ',CINIFILEPGD
718      WRITE(NLUOUT,FMT=*) ' whereas in NAM_REAL_PGD you have CPGD_FILE = '&
719           ,CPGD_FILE
720      WRITE(NLUOUT,FMT=*) ' '
721      WRITE(NLUOUT,FMT=*) ' CINIFILEPGD HAS BEEN SET TO  ',CPGD_FILE        
722      CINIFILEPGD=CPGD_FILE
723   END IF
724   IF ( IJPHEXT .NE. JPHEXT ) THEN
725      WRITE(NLUOUT,FMT=*) ' PREP_IDEAL_CASE : JPHEXT in PRE_IDEA1.nam/NAM_CONF_PRE ( or default value )&
726           JPHEXT=',JPHEXT
727      WRITE(NLUOUT,FMT=*) ' different from PGD files=', CINIFILEPGD,' value JPHEXT=',IJPHEXT
728      WRITE(NLUOUT,FMT=*) '-> JOB ABORTED'
729      CALL CLOSE_ll(CLUOUT,IOSTAT=IRESP)
730      CALL ABORT  
731      STOP   
732      !WRITE(NLUOUT,FMT=*) ' JPHEXT HAS BEEN SET TO ', IJPHEXT
733      !IJPHEXT = JPHEXT
734   END IF
735 END IF
736 !
737 NIMAX_ll=NIMAX   !! _ll variables are global variables
738 NJMAX_ll=NJMAX   !! but the old names are kept in PRE_IDEA1.nam file
739 !
740 !*       3.3   check some parameters:
741 !
742 L1D=.FALSE. ; L2D=.FALSE.
743 !
744 IF ((NIMAX == 1).OR.(NJMAX == 1)) THEN 
745   L2D=.TRUE.
746   NJMAX_ll=1
747   NIMAX_ll=MAX(NIMAX,NJMAX)
748   WRITE(NLUOUT,FMT=*) ' NJMAX HAS BEEN SET TO 1 SINCE 2D INITIAL FILE IS REQUIRED &
749                    & (L2D=TRUE) )' 
750 END IF
751 !
752 IF ((NIMAX == 1).AND.(NJMAX == 1)) THEN 
753   L1D=.TRUE.
754   NIMAX_ll = 1
755   NJMAX_ll = 1
756   WRITE(NLUOUT,FMT=*) ' 1D INITIAL FILE IS REQUIRED (L1D=TRUE) ' 
757 END IF
758 !
759 IF(.NOT. L1D) THEN
760   LHORELAX_UVWTH=.TRUE.
761   LHORELAX_RV=.TRUE.
762 ENDIF
763 !
764 NRIMX= MIN(JPRIMMAX,NIMAX_ll/2)
765 !
766 IF (L2D) THEN
767   NRIMY=0
768 ELSE
769   NRIMY= MIN(JPRIMMAX,NJMAX_ll/2)
770 END IF
771 !
772 IF (L1D) THEN
773   NRIMX=0
774   NRIMY=0
775 END IF
776 !
777 IF (L1D .AND. ( LPERTURB .OR. LGEOSBAL .OR.                &
778                (.NOT. LCARTESIAN ) .OR. (.NOT. LTHINSHELL) ))THEN 
779   LGEOSBAL   = .FALSE.
780   LPERTURB   = .FALSE.
781   LCARTESIAN = .TRUE. 
782   LTHINSHELL = .TRUE. 
783   WRITE(NLUOUT,FMT=*) ' LGEOSBAL AND LPERTURB HAVE BEEN SET TO FALSE &
784                       & AND LCARTESIAN AND LTHINSHELL TO TRUE        &
785                       & SINCE 1D INITIAL FILE IS REQUIRED (L1D=TRUE)' 
786 END IF
787 !
788 IF (LGEOSBAL .AND. LSHIFT ) THEN
789   LSHIFT=.FALSE.
790   WRITE(NLUOUT,FMT=*) ' LSHIFT HAS BEEN SET TO FALSE SINCE &
791                         & LGEOSBAL=.TRUE. IS REQUIRED '
792 END IF
793 !                      
794 !*       3.4   compute the number of moist variables :
795 !
796 IF (.NOT.LUSERV) THEN
797   LUSERV = .TRUE.
798   WRITE(NLUOUT,FMT=*) ' LUSERV HAS BEEN RESET TO TRUE, SINCE A MOIST VARIABLE &
799                    & IS PRESENT IN EXPRE FILE (CIDEAL = RSOU OR CSTN)' 
800 END IF
801 !
802 IF((LUSERI .OR. LUSERC).AND. (CIDEAL /= 'RSOU')) THEN
803   WRITE(NLUOUT,FMT=*) 'USE OF HYDROMETEORS IS ONLY ALLOWED IN RSOU CASE'
804   WRITE(NLUOUT,FMT=*) '-> JOB ABORTED'
805    !callabortstop
806    CALL CLOSE_ll(CLUOUT,IOSTAT=IRESP)
807    CALL ABORT
808   STOP
809 ENDIF
810 IF (LUSERI) THEN
811   LUSERC =.TRUE.
812   LUSERR =.TRUE.
813   LUSERI =.TRUE.
814   LUSERS =.TRUE.
815   LUSERG =.TRUE.
816   LUSERH =.FALSE.
817   CCLOUD='ICE3'
818 ELSEIF(LUSERC) THEN
819   LUSERR =.FALSE.
820   LUSERI =.FALSE.
821   LUSERS =.FALSE.
822   LUSERG =.FALSE.
823   LUSERH =.FALSE.
824   CCLOUD='REVE'
825 ELSE
826   LUSERC =.FALSE.
827   LUSERR =.FALSE.
828   LUSERI =.FALSE.
829   LUSERS =.FALSE.
830   LUSERG =.FALSE.
831   LUSERH =.FALSE.
832   LHORELAX_RC=.FALSE.
833   LHORELAX_RR=.FALSE.
834   LHORELAX_RI=.FALSE.
835   LHORELAX_RS=.FALSE.
836   LHORELAX_RG=.FALSE.
837   LHORELAX_RH=.FALSE.
838   CCLOUD='NONE'
839 !
840 END IF
841 !
842 NRR=0
843 IF (LUSERV) NRR=NRR+1
844 IF (LUSERC) NRR=NRR+1
845 IF (LUSERR) NRR=NRR+1
846 IF (LUSERI) NRR=NRR+1
847 IF (LUSERS) NRR=NRR+1
848 IF (LUSERG) NRR=NRR+1
849 IF (LUSERH) NRR=NRR+1
850 !
851 ! NRR=4 for RSOU case because RI and Rc always computed
852 IF (CIDEAL == 'RSOU' .AND. NRR < 4 ) NRR=4
853 !                   
854 !
855 !*       3.5   Chemistry
856 !
857 IF (LORILAM .OR. LCH_INIT_FIELD) THEN
858   ! Always initialize chemical scheme variables before INI_NSV call !
859   CALL CH_INIT_SCHEME_n(1,LUSECHAQ,LUSECHIC,LCH_PH,NLUOUT,NVERB)
860   LUSECHEM = .TRUE.
861   IF (LORILAM) THEN
862     CORGANIC = "MPMPO"
863     LVARSIGI = .TRUE.
864     LVARSIGJ = .TRUE.
865     CALL CH_AER_INIT_SOA(NLUOUT, NVERB)
866   END IF
867 END IF
868 ! initialise NSV_* variables
869 CALL INI_NSV(1)
870 LHORELAX_SV(:)=.FALSE.
871 IF(.NOT. L1D) LHORELAX_SV(1:NSV)=.TRUE.
872 !
873 !-------------------------------------------------------------------------------
874 !
875 !*       4.    ALLOCATE MEMORY FOR ARRAYS :  
876 !              ----------------------------
877 !
878 !*       4.1  Vertical Spatial grid 
879 !
880 CALL READ_VER_GRID(CEXPRE)
881 !
882 !*       4.2  Initialize parallel variables and compute array's dimensions
883 !
884 !
885 IF(LGEOSBAL) THEN
886   CALL SET_SPLITTING_ll('XSPLITTING')  ! required for integration of thermal wind balance
887 ELSE
888   CALL SET_SPLITTING_ll('BSPLITTING')
889 ENDIF
890 CALL SET_JP_ll(1,JPHEXT,JPVEXT,JPHEXT)
891 CALL SET_DAD0_ll()
892 CALL SET_DIM_ll(NIMAX_ll, NJMAX_ll, NKMAX)
893 CALL SET_FMPACK_ll(L1D,L2D,LPACK)
894 CALL SET_LBX_ll(CLBCX(1), 1)
895 CALL SET_LBY_ll(CLBCY(1), 1)
896 CALL SET_XRATIO_ll(1, 1)
897 CALL SET_YRATIO_ll(1, 1)
898 CALL SET_XOR_ll(1, 1)
899 CALL SET_XEND_ll(NIMAX_ll+2*JPHEXT, 1)
900 CALL SET_YOR_ll(1, 1)
901 CALL SET_YEND_ll(NJMAX_ll+2*JPHEXT, 1)
902 CALL SET_DAD_ll(0, 1)
903 CALL INI_PARAZ_ll(IINFO_ll)
904 !
905 ! sizes of arrays of the extended sub-domain
906 !
907 CALL GET_DIM_EXT_ll('B',NIU,NJU)
908 CALL GET_DIM_PHYS_ll('B',NIMAX,NJMAX)
909 CALL GET_INDICE_ll(NIB,NJB,NIE,NJE)
910 CALL GET_OR_ll('B',IXOR,IYOR)
911 NKB=1+JPVEXT
912 NKU=NKMAX+2*JPVEXT
913 !
914 !*       4.3  Global variables absent from the modules :
915 !
916 ALLOCATE(XJ(NIU,NJU,NKU))
917 SELECT CASE(CIDEAL)
918   CASE('RSOU','CSTN')
919     IF (LGEOSBAL) ALLOCATE(XCORIOZ(NIU,NJU,NKU))  ! exceptionally a 3D array  
920   CASE DEFAULT                      ! undefined preinitialization
921     WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE : CIDEAL IS NOT CORRECTLY DEFINED'
922     WRITE(NLUOUT,FMT=*) '-> JOB ABORTED'
923    !callabortstop
924    CALL CLOSE_ll(CLUOUT,IOSTAT=IRESP)
925    CALL ABORT
926     STOP
927 END SELECT 
928 !
929 !*       4.4   Prognostic variables at M instant (module MODD_FIELD1):
930 !
931 ALLOCATE(XUT(NIU,NJU,NKU))
932 ALLOCATE(XVT(NIU,NJU,NKU))
933 ALLOCATE(XWT(NIU,NJU,NKU))
934 ALLOCATE(XTHT(NIU,NJU,NKU))
935 ALLOCATE(XPABST(NIU,NJU,NKU))
936 ALLOCATE(XRT(NIU,NJU,NKU,NRR))
937 ALLOCATE(XSVT(NIU,NJU,NKU,NSV))
938 !
939 !*       4.5   Grid variables (module MODD_GRID1 and MODD_METRICS1):
940 !
941 ALLOCATE(XMAP(NIU,NJU))
942 ALLOCATE(XLAT(NIU,NJU))
943 ALLOCATE(XLON(NIU,NJU))
944 ALLOCATE(XDXHAT(NIU),XDYHAT(NJU))
945 IF (LEN_TRIM(CPGD_FILE)==0) ALLOCATE(XZS(NIU,NJU))
946 IF (LEN_TRIM(CPGD_FILE)==0) ALLOCATE(ZZS_ll(NIMAX_ll))
947 IF (LEN_TRIM(CPGD_FILE)==0) ALLOCATE(XZSMT(NIU,NJU))
948 ALLOCATE(XZZ(NIU,NJU,NKU))
949 !
950 ALLOCATE(XDXX(NIU,NJU,NKU))
951 ALLOCATE(XDYY(NIU,NJU,NKU))
952 ALLOCATE(XDZX(NIU,NJU,NKU))
953 ALLOCATE(XDZY(NIU,NJU,NKU))
954 ALLOCATE(XDZZ(NIU,NJU,NKU))
955 !
956 !*       4.6   Reference state variables (modules MODD_REF and MODD_REF1):
957 !
958 ALLOCATE(XRHODREFZ(NKU),XTHVREFZ(NKU))
959 XTHVREFZ(:)=0.0
960 IF(CEQNSYS == 'DUR') THEN
961   ALLOCATE(XRVREF(NIU,NJU,NKU))
962 ELSE
963   ALLOCATE(XRVREF(0,0,0))
964 END IF
965 ALLOCATE(XRHODREF(NIU,NJU,NKU),XTHVREF(NIU,NJU,NKU),XEXNREF(NIU,NJU,NKU))
966 ALLOCATE(XRHODJ(NIU,NJU,NKU))
967 !
968 !*       4.7   Larger Scale fields (modules MODD_LSFIELD1):
969 !
970 ALLOCATE(XLSUM(NIU,NJU,NKU))
971 ALLOCATE(XLSVM(NIU,NJU,NKU))
972 ALLOCATE(XLSWM(NIU,NJU,NKU))
973 ALLOCATE(XLSTHM(NIU,NJU,NKU))
974 ALLOCATE(XLSRVM(NIU,NJU,NKU))
975 !
976 !  allocate lateral boundary field used for coupling
977 !
978 IF ( L1D) THEN                         ! 1D case
979 !
980   NSIZELBX_ll=0
981   NSIZELBXU_ll=0
982   NSIZELBY_ll=0
983   NSIZELBYV_ll=0
984   NSIZELBXTKE_ll=0
985   NSIZELBXR_ll=0
986   NSIZELBXSV_ll=0
987   NSIZELBYTKE_ll=0
988   NSIZELBYR_ll=0
989   NSIZELBYSV_ll=0
990   ALLOCATE(XLBXUM(0,0,0))
991   ALLOCATE(XLBYUM(0,0,0))
992   ALLOCATE(XLBXVM(0,0,0))
993   ALLOCATE(XLBYVM(0,0,0))
994   ALLOCATE(XLBXWM(0,0,0))
995   ALLOCATE(XLBYWM(0,0,0))
996   ALLOCATE(XLBXTHM(0,0,0))
997   ALLOCATE(XLBYTHM(0,0,0))
998   ALLOCATE(XLBXTKEM(0,0,0))
999   ALLOCATE(XLBYTKEM(0,0,0))
1000   ALLOCATE(XLBXRM(0,0,0,0))
1001   ALLOCATE(XLBYRM(0,0,0,0))
1002   ALLOCATE(XLBXSVM(0,0,0,0))
1003   ALLOCATE(XLBYSVM(0,0,0,0))
1004 !
1005 ELSEIF( L2D ) THEN             ! 2D case (not yet parallelized)
1006 !                                          
1007   CALL GET_SIZEX_LB(CLUOUT,NIMAX_ll,NJMAX_ll,NRIMX,   &
1008        IISIZEXF,IJSIZEXF,IISIZEXFU,IJSIZEXFU,         &
1009        IISIZEX4,IJSIZEX4,IISIZEX2,IJSIZEX2)
1010   NSIZELBY_ll=0
1011   NSIZELBYV_ll=0
1012   NSIZELBYTKE_ll=0
1013   NSIZELBYR_ll=0
1014   NSIZELBYSV_ll=0
1015   ALLOCATE(XLBYUM(0,0,0))
1016   ALLOCATE(XLBYVM(0,0,0))
1017   ALLOCATE(XLBYWM(0,0,0))
1018   ALLOCATE(XLBYTHM(0,0,0))
1019   ALLOCATE(XLBYTKEM(0,0,0))
1020   ALLOCATE(XLBYRM(0,0,0,0))
1021   ALLOCATE(XLBYSVM(0,0,0,0))
1022   !
1023   IF ( LHORELAX_UVWTH ) THEN
1024 !JUAN A REVOIR TODO_JPHEXT
1025 ! <<<<<<< prep_ideal_case.f90
1026     ! NSIZELBX_ll=2*NRIMX+2
1027     ! NSIZELBXU_ll=2*NRIMX+2
1028     ALLOCATE(XLBXUM(IISIZEXFU,NJU,NKU))
1029     ALLOCATE(XLBXVM(IISIZEXF,NJU,NKU))
1030     ALLOCATE(XLBXWM(IISIZEXF,NJU,NKU))
1031     ALLOCATE(XLBXTHM(IISIZEXF,NJU,NKU))
1032 ! =======
1033     NSIZELBX_ll=2*NRIMX+2*JPHEXT
1034     NSIZELBXU_ll=2*NRIMX+2*JPHEXT
1035     ! ALLOCATE(XLBXUM(2*NRIMX+2*JPHEXT,NJU,NKU))
1036     ! ALLOCATE(XLBXVM(2*NRIMX+2*JPHEXT,NJU,NKU))
1037     ! ALLOCATE(XLBXWM(2*NRIMX+2*JPHEXT,NJU,NKU))
1038     ! ALLOCATE(XLBXTHM(2*NRIMX+2*JPHEXT,NJU,NKU))
1039 ! >>>>>>> 1.3.2.4.2.3.2.14.2.8.2.11.2.2
1040   ELSE
1041     NSIZELBX_ll= 2*JPHEXT     ! 2
1042     NSIZELBXU_ll=2*(JPHEXT+1) ! 4 
1043     ALLOCATE(XLBXUM(NSIZELBXU_ll,NJU,NKU))
1044     ALLOCATE(XLBXVM(NSIZELBX_ll,NJU,NKU))
1045     ALLOCATE(XLBXWM(NSIZELBX_ll,NJU,NKU))
1046     ALLOCATE(XLBXTHM(NSIZELBX_ll,NJU,NKU))
1047   END IF  
1048   !
1049   IF ( NRR > 0 ) THEN
1050     IF (       LHORELAX_RV .OR. LHORELAX_RC .OR. LHORELAX_RR .OR. LHORELAX_RI    &
1051           .OR. LHORELAX_RS .OR. LHORELAX_RG .OR. LHORELAX_RH                     &
1052        ) THEN 
1053 !JUAN A REVOIR TODO_JPHEXT
1054 ! <<<<<<< prep_ideal_case.f90
1055       ! NSIZELBXR_ll=2* NRIMX+2
1056       ALLOCATE(XLBXRM(IISIZEXF,NJU,NKU,NRR))
1057 ! =======
1058       NSIZELBXR_ll=2*NRIMX+2*JPHEXT
1059       ! ALLOCATE(XLBXRM(2*NRIMX+2*JPHEXT,NJU,NKU,NRR))
1060 ! >>>>>>> 1.3.2.4.2.3.2.14.2.8.2.11.2.2
1061     ELSE
1062       NSIZELBXR_ll=2*JPHEXT ! 2
1063       ALLOCATE(XLBXRM(NSIZELBXR_ll,NJU,NKU,NRR))
1064     ENDIF
1065   ELSE
1066     NSIZELBXR_ll=0
1067     ALLOCATE(XLBXRM(0,0,0,0))
1068   END IF
1069   !
1070   IF ( NSV > 0 ) THEN 
1071     IF ( ANY( LHORELAX_SV(:)) ) THEN
1072 !JUAN A REVOIR TODO_JPHEXT
1073 ! <<<<<<< prep_ideal_case.f90
1074       ! NSIZELBXSV_ll=2* NRIMX+2
1075       ALLOCATE(XLBXSVM(IISIZEXF,NJU,NKU,NSV))
1076 ! =======
1077       NSIZELBXSV_ll=2*NRIMX+2*JPHEXT
1078       ! ALLOCATE(XLBXSVM(2*NRIMX+2*JPHEXT,NJU,NKU,NSV))
1079 ! >>>>>>> 1.3.2.4.2.3.2.14.2.8.2.11.2.2
1080     ELSE
1081       NSIZELBXSV_ll=2*JPHEXT ! 2
1082       ALLOCATE(XLBXSVM(NSIZELBXSV_ll,NJU,NKU,NSV))
1083     END IF
1084   ELSE
1085     NSIZELBXSV_ll=0
1086     ALLOCATE(XLBXSVM(0,0,0,0))
1087   END IF
1088 !
1089 ELSE                                   ! 3D case
1090 !
1091   CALL GET_SIZEX_LB(CLUOUT,NIMAX_ll,NJMAX_ll,NRIMX,   &
1092        IISIZEXF,IJSIZEXF,IISIZEXFU,IJSIZEXFU,         &
1093        IISIZEX4,IJSIZEX4,IISIZEX2,IJSIZEX2)
1094   CALL GET_SIZEY_LB(CLUOUT,NIMAX_ll,NJMAX_ll,NRIMY,   &
1095        IISIZEYF,IJSIZEYF,IISIZEYFV,IJSIZEYFV,         &
1096        IISIZEY4,IJSIZEY4,IISIZEY2,IJSIZEY2)
1097 !
1098   IF ( LHORELAX_UVWTH ) THEN
1099     NSIZELBX_ll=2*NRIMX+2*JPHEXT
1100     NSIZELBXU_ll=2*NRIMX+2*JPHEXT
1101     NSIZELBY_ll=2*NRIMY+2*JPHEXT
1102     NSIZELBYV_ll=2*NRIMY+2*JPHEXT
1103     ALLOCATE(XLBXUM(IISIZEXFU,IJSIZEXFU,NKU))
1104     ALLOCATE(XLBYUM(IISIZEYF,IJSIZEYF,NKU))
1105     ALLOCATE(XLBXVM(IISIZEXF,IJSIZEXF,NKU))
1106     ALLOCATE(XLBYVM(IISIZEYFV,IJSIZEYFV,NKU))
1107     ALLOCATE(XLBXWM(IISIZEXF,IJSIZEXF,NKU))
1108     ALLOCATE(XLBYWM(IISIZEYF,IJSIZEYF,NKU))
1109     ALLOCATE(XLBXTHM(IISIZEXF,IJSIZEXF,NKU))
1110     ALLOCATE(XLBYTHM(IISIZEYF,IJSIZEYF,NKU))
1111   ELSE
1112     NSIZELBX_ll=2*JPHEXT      ! 2
1113     NSIZELBXU_ll=2*(JPHEXT+1) ! 4
1114     NSIZELBY_ll=2*JPHEXT      ! 2
1115     NSIZELBYV_ll=2*(JPHEXT+1) ! 4
1116     ALLOCATE(XLBXUM(IISIZEX4,IJSIZEX4,NKU))
1117     ALLOCATE(XLBYUM(IISIZEY2,IJSIZEY2,NKU))
1118     ALLOCATE(XLBXVM(IISIZEX2,IJSIZEX2,NKU))
1119     ALLOCATE(XLBYVM(IISIZEY4,IJSIZEY4,NKU))
1120     ALLOCATE(XLBXWM(IISIZEX2,IJSIZEX2,NKU))
1121     ALLOCATE(XLBYWM(IISIZEY2,IJSIZEY2,NKU))
1122     ALLOCATE(XLBXTHM(IISIZEX2,IJSIZEX2,NKU))
1123     ALLOCATE(XLBYTHM(IISIZEY2,IJSIZEY2,NKU))
1124   END IF  
1125   !
1126   IF ( NRR > 0 ) THEN
1127     IF (       LHORELAX_RV .OR. LHORELAX_RC .OR. LHORELAX_RR .OR. LHORELAX_RI    &
1128           .OR. LHORELAX_RS .OR. LHORELAX_RG .OR. LHORELAX_RH                     &
1129        ) THEN 
1130       NSIZELBXR_ll=2*NRIMX+2*JPHEXT
1131       NSIZELBYR_ll=2*NRIMY+2*JPHEXT
1132       ALLOCATE(XLBXRM(IISIZEXF,IJSIZEXF,NKU,NRR))
1133       ALLOCATE(XLBYRM(IISIZEYF,IJSIZEYF,NKU,NRR))
1134     ELSE
1135       NSIZELBXR_ll=2*JPHEXT    ! 2
1136       NSIZELBYR_ll=2*JPHEXT    ! 2
1137       ALLOCATE(XLBXRM(IISIZEX2,IJSIZEX2,NKU,NRR))
1138       ALLOCATE(XLBYRM(IISIZEY2,IJSIZEY2,NKU,NRR))
1139     ENDIF
1140   ELSE
1141     NSIZELBXR_ll=0
1142     NSIZELBYR_ll=0
1143     ALLOCATE(XLBXRM(0,0,0,0))
1144     ALLOCATE(XLBYRM(0,0,0,0))
1145   END IF
1146   !
1147   IF ( NSV > 0 ) THEN 
1148     IF ( ANY( LHORELAX_SV(:)) ) THEN
1149       NSIZELBXSV_ll=2*NRIMX+2*JPHEXT
1150       NSIZELBYSV_ll=2*NRIMY+2*JPHEXT
1151       ALLOCATE(XLBXSVM(IISIZEXF,IJSIZEXF,NKU,NSV))
1152       ALLOCATE(XLBYSVM(IISIZEYF,IJSIZEYF,NKU,NSV))
1153     ELSE
1154       NSIZELBXSV_ll=2*JPHEXT    ! 2
1155       NSIZELBYSV_ll=2*JPHEXT    ! 2
1156       ALLOCATE(XLBXSVM(IISIZEX2,IJSIZEX2,NKU,NSV))
1157       ALLOCATE(XLBYSVM(IISIZEY2,IJSIZEY2,NKU,NSV))
1158     END IF
1159   ELSE
1160     NSIZELBXSV_ll=0
1161     NSIZELBYSV_ll=0
1162     ALLOCATE(XLBXSVM(0,0,0,0))
1163     ALLOCATE(XLBYSVM(0,0,0,0))
1164   END IF
1165 END IF
1166 !
1167 !
1168 !-------------------------------------------------------------------------------
1169 !
1170 !*       5.     INITIALIZE ALL THE MODEL VARIABLES
1171 !               ----------------------------------
1172 !
1173 !
1174 !*       5.1    Grid variables and RS localization:
1175 !
1176 !*       5.1.1  Horizontal Spatial grid :
1177 !
1178 IF( LEN_TRIM(CPGD_FILE) /= 0 ) THEN 
1179 !--------------------------------------------------------
1180 ! the MESONH horizontal grid will be read in the PGD_FILE 
1181 !--------------------------------------------------------
1182   CALL READ_HGRID(1,CPGD_FILE,YPGD_NAME,YPGD_DAD_NAME,YPGD_TYPE)
1183 ! control the cartesian option
1184   IF( LCARTESIAN ) THEN
1185      WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE : IN GENERAL, THE USE OF A PGD_FILE &
1186                 & IMPLIES THAT YOU MUST TAKE INTO ACCOUNT THE EARTH SPHERICITY'
1187      WRITE(NLUOUT,FMT=*) 'NEVERTHELESS, LCARTESIAN HAS BEEN KEPT TO TRUE'
1188   END IF   
1189 !
1190 !* use of the externalized surface
1191 !
1192   CSURF = "EXTE"
1193 !
1194 ! determine whether the model is flat or no
1195 !
1196   ZZS_MAX = ABS( MAXVAL(XZS(NIB:NIU-JPHEXT,NJB:NJU-JPHEXT)))
1197   CALL MPI_ALLREDUCE(ZZS_MAX, ZZS_MAX_ll, 1, MPI_PRECISION, MPI_MAX,  &
1198                      NMNH_COMM_WORLD,IINFO_ll)
1199   IF( ABS(ZZS_MAX_ll)  < 1.E-10 ) THEN
1200     LFLAT=.TRUE.
1201   ELSE
1202     LFLAT=.FALSE.
1203   END IF
1204 !
1205
1206 ELSE
1207 !------------------------------------------------------------------------
1208 ! the MESONH horizontal grid is built from the PRE_IDEA1.nam informations
1209 !------------------------------------------------------------------------
1210 !
1211   ALLOCATE(XXHAT(NIU),XYHAT(NJU))
1212 !
1213 ! define the grid localization at the earth surface by the central point
1214 ! coordinates
1215 !
1216   IF (XLONCEN/=XUNDEF .OR. XLATCEN/=XUNDEF) THEN
1217     IF (XLONCEN/=XUNDEF .AND. XLATCEN/=XUNDEF) THEN 
1218 !
1219 ! it should be noted that XLATCEN and XLONCEN refer to a vertical
1220 ! vorticity point and (XLATORI, XLONORI) refer to the mass point of
1221 ! conformal coordinates (0,0). This is to allow the centering of the model in
1222 ! a non-cyclic  configuration regarding to XLATCEN or XLONCEN.
1223 !
1224       ALLOCATE(ZXHAT_ll(NIMAX_ll+2*JPHEXT),ZYHAT_ll(NJMAX_ll+2*JPHEXT))
1225       ZXHAT_ll=0.
1226       ZYHAT_ll=0.
1227       CALL SM_LATLON(XLATCEN,XLONCEN,                     &
1228                        -XDELTAX*(NIMAX_ll/2-0.5+JPHEXT),  &
1229                        -XDELTAY*(NJMAX_ll/2-0.5+JPHEXT),  &
1230                        XLATORI,XLONORI)
1231         DEALLOCATE(ZXHAT_ll,ZYHAT_ll)
1232 !
1233       WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE : XLATORI=' , XLATORI, &
1234                           ' XLONORI= ', XLONORI
1235     ELSE
1236       WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE : LATITUDE AND LONGITUDE OF THE CENTER &
1237                            & POINT MUST BE INITIALIZED ALL TOGETHER OR NOT'
1238       WRITE(NLUOUT,FMT=*) '-> JOB ABORTED'
1239    !callabortstop
1240       CALL CLOSE_ll(CLUOUT,IOSTAT=IRESP)
1241       CALL ABORT
1242       STOP
1243     END IF
1244   END IF
1245 !
1246   IF (NPROC > 1) THEN
1247     CALL GET_DIM_EXT_ll('B',IXDIM,IYDIM)
1248     IBEG = IXOR-JPHEXT-1
1249     IEND = IBEG+IXDIM-1
1250     XXHAT(:) = (/ (FLOAT(JLOOP)*XDELTAX, JLOOP=IBEG,IEND) /)
1251     IBEG = IYOR-JPHEXT-1
1252     IEND = IBEG+IYDIM-1
1253     XYHAT(:) = (/ (FLOAT(JLOOP)*XDELTAY, JLOOP=IBEG,IEND) /)
1254 !
1255   ELSE
1256     XXHAT(:) = (/ (FLOAT(JLOOP-NIB)*XDELTAX, JLOOP=1,NIU) /)
1257     XYHAT(:) = (/ (FLOAT(JLOOP-NJB)*XDELTAY, JLOOP=1,NJU) /)
1258   END IF
1259 END IF
1260 !
1261 !*       5.1.2  Orography and Gal-Chen Sommerville transformation :
1262 !
1263 IF (    LEN_TRIM(CPGD_FILE) == 0  .OR. .NOT. LREAD_ZS) THEN
1264   SELECT CASE(CZS)     ! 'FLAT' or 'SINE' or 'BELL'
1265   CASE('FLAT')
1266     LFLAT = .TRUE.
1267     IF (XHMAX==XUNDEF) THEN
1268       XZS(:,:) = 0.
1269     ELSE
1270       XZS(:,:) = XHMAX
1271     END IF
1272   CASE('SINE')       ! sinus-shaped orography 
1273     IF (XHMAX==XUNDEF) XHMAX=300.
1274     LFLAT    =.FALSE.
1275     XZS(:,:) = XHMAX          &      ! three-dimensional case   
1276     *SPREAD((/((SIN((XPI/(NIMAX_ll+2*JPHEXT-1))*JLOOP)**2)**NEXPX,JLOOP=IXOR-1,IXOR+NIU-2)/),2,NJU) &
1277     *SPREAD((/((SIN((XPI/(NJMAX_ll+2*JPHEXT-1))*JLOOP)**2)**NEXPY,JLOOP=IYOR-1,IYOR+NJU-2)/),1,NIU)
1278     IF(L1D) THEN                     ! one-dimensional case
1279       XZS(:,:) = XHMAX 
1280     END IF        
1281   CASE('BELL')       ! bell-shaped orography 
1282     IF (XHMAX==XUNDEF) XHMAX=300.
1283     LFLAT = .FALSE.
1284     IF(.NOT.L2D) THEN                ! three-dimensional case
1285       XZS(:,:) = XHMAX  / ( 1.                                           &
1286         + ( (SPREAD(XXHAT(1:NIU),2,NJU) - FLOAT(NIZS) * XDELTAX) /XAX ) **2  &
1287         + ( (SPREAD(XYHAT(1:NJU),1,NIU) - FLOAT(NJZS) * XDELTAY) /XAY ) **2  ) **1.5
1288     ELSE                             ! two-dimensional case
1289       XZS(:,:) = XHMAX  / ( 1.                                          &
1290         + ( (SPREAD(XXHAT(1:NIU),2,NJU) - FLOAT(NIZS) * XDELTAX) /XAX ) **2 )
1291     ENDIF
1292     IF(L1D) THEN                     ! one-dimensional case
1293       XZS(:,:) = XHMAX 
1294     END IF        
1295   CASE('COSI')       ! (1+cosine)**4 shape
1296     IF (XHMAX==XUNDEF) XHMAX=800.
1297     LFLAT = .FALSE.
1298     IF(L2D) THEN                     ! two-dimensional case
1299       DO JILOOP = 1, NIU
1300         ZDIST = XXHAT(JILOOP)-FLOAT(NIZS)*XDELTAX
1301         IF( ABS(ZDIST)<(4.0*XAX) ) THEN
1302           XZS(JILOOP,:) = (XHMAX/16.0)*( 1.0 + COS((XPI*ZDIST)/(4.0*XAX)) )**4
1303         ELSE
1304           XZS(JILOOP,:) = 0.0
1305         ENDIF
1306       END DO
1307     ENDIF
1308   CASE('SCHA')       ! exp(-(x/a)**2)*cosine(pi*x/lambda)**2 shape
1309     IF (XHMAX==XUNDEF) XHMAX=800.
1310     LFLAT = .FALSE.
1311     IF(L2D) THEN                     ! two-dimensional case
1312       DO JILOOP = 1, NIU
1313         ZDIST = XXHAT(JILOOP)-FLOAT(NIZS)*XDELTAX
1314         IF( ABS(ZDIST)<(4.0*XAX) ) THEN
1315           XZS(JILOOP,:) = XHMAX*EXP(-(ZDIST/XAY)**2)*COS((XPI*ZDIST)/XAX)**2
1316         ELSE
1317           XZS(JILOOP,:) = 0.0
1318         ENDIF
1319       END DO
1320     ENDIF
1321   CASE('AGNE')       ! h*a**2/(x**2+a**2) shape
1322     LFLAT = .FALSE.
1323     IF(L2D) THEN                     ! two-dimensional case
1324       DO JILOOP = 1, NIU
1325         ZDIST = XXHAT(JILOOP)-FLOAT(NIZS)*XDELTAX
1326           XZS(JILOOP,:) = XHMAX*(XAX**2)/(XAX**2+ZDIST**2)
1327       END DO
1328                 ELSE            ! three dimensionnal case - infinite profile in y direction
1329                         DO JILOOP = 1, NIU
1330         ZDIST = XXHAT(JILOOP)-FLOAT(NIZS)*XDELTAX
1331           XZS(JILOOP,:) = XHMAX*(XAX**2)/(XAX**2+ZDIST**2)
1332       END DO
1333     ENDIF
1334
1335   CASE('DATA')       ! discretized orography
1336     LFLAT    =.FALSE.
1337     WRITE(NLUOUT,FMT=*) 'CZS="DATA",   ATTEMPT TO READ ARRAY     &
1338                     &XZS(NIB:NIU-JPHEXT:1,NJU-JPHEXT:NJB:-1) &
1339                     &starting from the first index'
1340     CALL POSKEY(NLUPRE,NLUOUT,'ZSDATA')
1341     DO JJLOOP = NJMAX_ll+2*JPHEXT-1,JPHEXT+1,-1    ! input like a map prior the sounding
1342       READ(NLUPRE,FMT=*) ZZS_ll
1343       IF ( ( JJLOOP <= ( NJU-JPHEXT + IYOR-1 ) ) .AND. ( JJLOOP >= ( NJB + IYOR-1 ) ) ) THEN
1344          IJ    = JJLOOP - ( IYOR-1 )
1345          XZS(NIB:NIU-JPHEXT,IJ) = ZZS_ll(IXOR:IXOR + NIU-JPHEXT - NIB )
1346       END IF
1347     END DO
1348 !
1349   CASE DEFAULT   ! undefined  shape of orography
1350     WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE: ERRONEOUS TERRAIN TYPE'
1351     WRITE(NLUOUT,FMT=*) '-> JOB ABORTED'
1352    !callabortstop
1353     CALL CLOSE_ll(CLUOUT,IOSTAT=IRESP)
1354     CALL ABORT
1355     STOP
1356   END SELECT
1357 !
1358   CALL ADD2DFIELD_ll(TZ_FIELDS_ll, XZS)
1359   CALL UPDATE_HALO_ll(TZ_FIELDS_ll,IINFO_ll)
1360   CALL CLEANLIST_ll(TZ_FIELDS_ll)
1361 !
1362 END IF
1363 !
1364 !IF( ( LEN_TRIM(CPGD_FILE) /= 0 ) .AND. .NOT.LFLAT .AND. &
1365 ! ((CLBCX(1) /= "OPEN" ) .OR. &
1366 ! (CLBCX(2) /= "OPEN" ) .OR. (CLBCY(1) /= "OPEN" ) .OR. &
1367 ! (CLBCY(2) /= "OPEN" )) )  THEN 
1368 !  WRITE(NLUOUT,FMT=*) 'STOP:WITH A PGD FILE YOU CANNOT BE IN CYCLIC LBC'
1369 !   !callabortstop
1370 !  CALL CLOSE_ll(CLUOUT,IOSTAT=IRESP)
1371 !  CALL ABORT
1372 !  STOP
1373 !END IF
1374 !
1375 IF (LWEST_ll())  THEN
1376   DO JILOOP = 1,JPHEXT
1377     XZS(JILOOP,:) = XZS(NIB,:)
1378   END DO
1379 END IF
1380 IF (LEAST_ll()) THEN
1381   DO JILOOP = NIU-JPHEXT+1,NIU
1382     XZS(JILOOP,:)=XZS(NIU-JPHEXT,:)
1383   END DO
1384 END IF
1385 IF (LSOUTH_ll()) THEN
1386   DO JJLOOP = 1,JPHEXT
1387     XZS(:,JJLOOP)=XZS(:,NJB)
1388   END DO
1389 END IF
1390 IF (LNORTH_ll()) THEN
1391   DO JJLOOP =NJU-JPHEXT+1,NJU
1392     XZS(:,JJLOOP)=XZS(:,NJU-JPHEXT)
1393   END DO
1394 END IF
1395 !
1396 IF ( LEN_TRIM(CPGD_FILE) == 0  .OR. .NOT. LREAD_ZS) THEN
1397   IF (LSLEVE) THEN
1398     CALL ZSMT_PIC(NSLEVE,XSMOOTH_ZS)
1399   ELSE
1400     XZSMT(:,:) = 0.
1401   END IF
1402 END IF
1403 !
1404 IF (LCARTESIAN) THEN
1405   CALL SM_GRIDCART(CLUOUT,XXHAT,XYHAT,XZHAT,XZS,LSLEVE,XLEN1,XLEN2,XZSMT,XDXHAT,XDYHAT,XZZ,XJ)
1406   XMAP=1.
1407 ELSE
1408   CALL SM_GRIDPROJ(CLUOUT,XXHAT,XYHAT,XZHAT,XZS,LSLEVE,XLEN1,XLEN2,XZSMT,XLATORI,XLONORI, &
1409                    XMAP,XLAT,XLON,XDXHAT,XDYHAT,XZZ,XJ)
1410 END IF
1411 !*       5.4.1  metrics coefficients and update halos:
1412 !
1413 CALL METRICS(XMAP,XDXHAT,XDYHAT,XZZ,XDXX,XDYY,XDZX,XDZY,XDZZ)
1414 !
1415 CALL UPDATE_METRICS(CLBCX,CLBCY,XDXX,XDYY,XDZX,XDZY,XDZZ)                         
1416 !
1417 !*       5.1.3  Compute the localization in index space of the vertical profile
1418 !               in CSTN and RSOU cases  :
1419 !
1420 IF (CTYPELOC =='LATLON' ) THEN  
1421   IF (.NOT.LCARTESIAN) THEN                            ! compute (x,y) if 
1422     CALL SM_XYHAT(XLATORI,XLONORI,                 &   ! the localization 
1423                   XLATLOC,XLONLOC,XXHATLOC,XYHATLOC)   ! is given in latitude 
1424   ELSE                                                 ! and longitude
1425     WRITE(NLUOUT,FMT=*) 'CTYPELOC CANNOT BE LATLON IN CARTESIAN GEOMETRY'
1426     WRITE(NLUOUT,FMT=*) '-> JOB ABORTED'
1427    !callabortstop
1428     CALL CLOSE_ll(CLUOUT,IOSTAT=IRESP)
1429     CALL ABORT
1430     STOP
1431   END IF 
1432 END IF  
1433 !
1434 ALLOCATE(ZXHAT_ll(NIMAX_ll+ 2 * JPHEXT),ZYHAT_ll(NJMAX_ll+2 * JPHEXT))
1435 CALL GATHERALL_FIELD_ll('XX',XXHAT,ZXHAT_ll,NRESP) !//
1436 CALL GATHERALL_FIELD_ll('YY',XYHAT,ZYHAT_ll,NRESP) !//
1437 IF (CTYPELOC /= 'IJGRID') THEN                                               
1438   NILOC = MINLOC(ABS(XXHATLOC-ZXHAT_ll(:)))
1439   NJLOC = MINLOC(ABS(XYHATLOC-ZYHAT_ll(:)))
1440 END IF
1441 !
1442 IF ( L1D .AND. ( NILOC(1) /= 1 .OR. NJLOC(1) /= 1 ) ) THEN
1443   NILOC = 1
1444   NJLOC = 1
1445   WRITE(NLUOUT,FMT=*) 'FOR 1D CONFIGURATION, THE RS INFORMATIONS ARE TAKEN AT &
1446                       & I=1 AND J=1 (CENTRAL VERTICAL WITHOUT HALO)'
1447 END IF
1448 !
1449 IF ( L2D .AND. ( NJLOC(1) /= 1 ) ) THEN
1450   NJLOC = 1
1451   WRITE(NLUOUT,FMT=*) 'FOR 2D CONFIGURATION, THE RS INFORMATIONS ARE TAKEN AT &
1452                       & J=1 (CENTRAL PLANE WITHOUT HALO)'
1453 END IF
1454 !
1455 !*       5.2    Prognostic variables (not multiplied by  rhoJ) : u,v,w,theta,r
1456 !               and 1D anelastic reference state
1457 !
1458 IF(LPV_PERT .AND. .NOT.(LGEOSBAL)) THEN
1459   WRITE(NLUOUT,FMT=*) 'FOR PV INVERSION, LGEOSBAL HAS TO BE TRUE'
1460    !callabortstop
1461   CALL CLOSE_ll(CLUOUT,IOSTAT=IRESP)
1462   CALL ABORT
1463   STOP
1464 ENDIF
1465 !
1466 IF(LPV_PERT .AND. NPROC>1) THEN
1467     WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE : THE USE OF A PV INVERSION HAS TO BE &
1468                         & PERFORMED WITH MONOPROCESSOR MODE'
1469     WRITE(NLUOUT,FMT=*) '-> JOB ABORTED'
1470    !callabortstop
1471    CALL CLOSE_ll(CLUOUT,IOSTAT=IRESP)
1472    CALL ABORT
1473     STOP
1474 ENDIF
1475 !
1476 !*       5.2.1  Use a Radiosounding : CIDEAL='RSOU''
1477 !
1478 IF (CIDEAL == 'RSOU') THEN
1479   WRITE(NLUOUT,FMT=*) 'CIDEAL="RSOU", attempt to read DATE'
1480   CALL POSKEY(NLUPRE,NLUOUT,'RSOU')
1481   READ(NLUPRE,FMT=*)  NYEAR,NMONTH,NDAY,XTIME
1482   TDTCUR = DATE_TIME(DATE(NYEAR,NMONTH,NDAY),XTIME)
1483   TDTEXP = TDTCUR
1484   TDTSEG = TDTCUR
1485   TDTMOD = TDTCUR
1486   READ(NLUPRE,*) YKIND
1487   BACKSPACE(NLUPRE)    ! because YKIND read again in set_rsou
1488   WRITE(NLUOUT,FMT=*) 'CIDEAL="RSOU", ATTEMPT TO PROCESS THE SOUNDING DATA'
1489   IF (LGEOSBAL) THEN
1490     CALL SET_RSOU(CEXPRE,CFUNU,CFUNV,NILOC(1),NJLOC(1),LBOUSS,LPV_PERT,&
1491                   LRMV_BL,XJ,LSHIFT,XCORIOZ)
1492   ELSE
1493     CALL SET_RSOU(CEXPRE,CFUNU,CFUNV,NILOC(1),NJLOC(1),LBOUSS,LPV_PERT,&
1494                   LRMV_BL,XJ,LSHIFT)
1495   END IF
1496 !
1497 !*       5.2.2  N=cste  and U(z) : CIDEAL='CSTN'
1498 !
1499 ELSE IF (CIDEAL == 'CSTN') THEN
1500   WRITE(NLUOUT,FMT=*) 'CIDEAL="CSTN", attempt to read DATE'
1501   CALL POSKEY(NLUPRE,NLUOUT,'CSTN')
1502   READ(NLUPRE,FMT=*)  NYEAR,NMONTH,NDAY,XTIME
1503   TDTCUR = DATE_TIME(DATE(NYEAR,NMONTH,NDAY),XTIME)
1504   TDTEXP = TDTCUR
1505   TDTSEG = TDTCUR
1506   TDTMOD = TDTCUR
1507   WRITE(NLUOUT,FMT=*) 'CIDEAL="CSTN", ATTEMPT TO PROCESS THE SOUNDING DATA'
1508   IF (LGEOSBAL) THEN
1509     CALL SET_CSTN(CEXPRE,CFUNU,CFUNV,NILOC(1),NJLOC(1),LBOUSS,LPV_PERT,&
1510                   LRMV_BL,XJ,LSHIFT,XCORIOZ)
1511   ELSE
1512     CALL SET_CSTN(CEXPRE,CFUNU,CFUNV,NILOC(1),NJLOC(1),LBOUSS,LPV_PERT,&
1513                   LRMV_BL,XJ,LSHIFT)
1514   END IF
1515 !
1516 END IF 
1517 !
1518 !*       5.3    Forcing variables
1519 !
1520 IF (LFORCING) THEN
1521   WRITE(NLUOUT,FMT=*) 'FORCING IS ENABLED, ATTEMPT TO SET FORCING FIELDS'
1522   CALL POSKEY(NLUPRE,NLUOUT,'ZFRC ','PFRC')
1523   CALL SET_FRC(CEXPRE)
1524 END IF
1525 !
1526 !! ---------------------------------------------------------------------
1527 ! Modif PP ADV FRC
1528 ! 5.4.2 initialize profiles for adv forcings
1529 IF (L2D_ADV_FRC) THEN
1530     WRITE(NLUOUT,FMT=*) 'L2D_ADV_FRC IS SET TO  TRUE'
1531     WRITE(NLUOUT,FMT=*) 'ADVECTING FORCING USED IS USER MADE, NOT STANDARD ONE ' 
1532     WRITE(NLUOUT,FMT=*) 'IT IS FOR 2D IDEALIZED WAM STUDY ONLY ' 
1533    CALL POSKEY(NLUPRE,NLUOUT,'ZFRC_ADV')
1534    CALL SET_ADVFRC(CEXPRE)
1535 ENDIF
1536 IF (L2D_REL_FRC) THEN
1537     WRITE(NLUOUT,FMT=*) 'L2D_REL_FRC IS SET TO  TRUE'
1538     WRITE(NLUOUT,FMT=*) 'RELAXATION FORCING USED IS USER MADE, NOT STANDARD ONE ' 
1539     WRITE(NLUOUT,FMT=*) 'IT IS FOR 2D IDEALIZED WAM STUDY ONLY ' 
1540    CALL POSKEY(NLUPRE,NLUOUT,'ZFRC_REL')
1541    CALL SET_RELFRC(CEXPRE)
1542 ENDIF
1543 !*       5.4    3D Reference state variables :
1544 !
1545 !
1546 !*       5.4.1  metrics coefficients and update halos:
1547 !
1548 CALL METRICS(XMAP,XDXHAT,XDYHAT,XZZ,XDXX,XDYY,XDZX,XDZY,XDZZ)
1549 !
1550 CALL UPDATE_METRICS(CLBCX,CLBCY,XDXX,XDYY,XDZX,XDZY,XDZZ)
1551 !
1552 !*       5.4.2  3D reference state :
1553 !
1554 CALL SET_REF(0,'NIL',CLUOUT,                                     &
1555              XZZ,XZHAT,XJ,XDXX,XDYY,CLBCX,CLBCY,                 &
1556              XREFMASS,XMASS_O_PHI0,XLINMASS,                     &
1557              XRHODREF,XTHVREF,XRVREF,XEXNREF,XRHODJ              )
1558 !
1559 !
1560 !*       5.5.1  Absolute pressure :
1561 !
1562 !
1563 !*       5.5.2  Total mass of dry air Md computation :
1564 !
1565 CALL TOTAL_DMASS(CLUOUT,XJ,XRHODREF,XDRYMASST)
1566 !
1567 !
1568 !*       5.6    Complete prognostic variables (multipliy by  rhoJ) at time t :
1569 !
1570 ! U grid   : gridpoint 2
1571 IF (LWEST_ll())  XUT(1,:,:)    = 2.*XUT(2,:,:) - XUT(3,:,:)
1572 ! V grid   : gridpoint 3
1573 IF (LSOUTH_ll())  XVT(:,1,:)    = 2.*XVT(:,2,:) - XVT(:,3,:)
1574 ! SV : gridpoint 1
1575 XSVT(:,:,:,:) = 0.
1576 !
1577 !
1578 !*       5.7   Larger scale fields initialization :
1579 !
1580 XLSUM(:,:,:) = XUT(:,:,:)        ! these fields do not satisfy the 
1581 XLSVM(:,:,:) = XVT(:,:,:)        ! lower boundary condition but are 
1582 XLSWM(:,:,:) = XWT(:,:,:)        ! in equilibrium
1583 XLSTHM(:,:,:)= XTHT(:,:,:)
1584 XLSRVM(:,:,:)= XRT(:,:,:,1)
1585 !
1586 ! enforce the vertical homogeneity under the ground and above the top of
1587 ! the model for the LS fields
1588 !
1589 XLSUM(:,:,NKB-1)=XLSUM(:,:,NKB)
1590 XLSUM(:,:,NKU)=XLSUM(:,:,NKU-1)
1591 XLSVM(:,:,NKB-1)=XLSVM(:,:,NKB)
1592 XLSVM(:,:,NKU)=XLSVM(:,:,NKU-1)
1593 XLSWM(:,:,NKB-1)=XLSWM(:,:,NKB)
1594 XLSWM(:,:,NKU)=XLSWM(:,:,NKU-1)
1595 XLSTHM(:,:,NKB-1)=XLSTHM(:,:,NKB)
1596 XLSTHM(:,:,NKU)=XLSTHM(:,:,NKU-1)
1597 IF ( NRR > 0 ) THEN
1598   XLSRVM(:,:,NKB-1)=XLSRVM(:,:,NKB)
1599   XLSRVM(:,:,NKU)=XLSRVM(:,:,NKU-1)
1600 END IF
1601 !
1602 ILBX=SIZE(XLBXUM,1)
1603 ILBY=SIZE(XLBYUM,2)
1604 IF(LWEST_ll() .AND. .NOT. L1D) THEN
1605   XLBXUM(1:NRIMX+JPHEXT,        :,:)     = XUT(2:NRIMX+JPHEXT+1,        :,:)
1606   XLBXVM(1:NRIMX+JPHEXT,        :,:)     = XVT(1:NRIMX+JPHEXT,        :,:)
1607   XLBXWM(1:NRIMX+JPHEXT,        :,:)     = XWT(1:NRIMX+JPHEXT,        :,:)
1608   XLBXTHM(1:NRIMX+JPHEXT,        :,:)   = XTHT(1:NRIMX+JPHEXT,        :,:)
1609   XLBXRM(1:NRIMX+JPHEXT,        :,:,:)   = XRT(1:NRIMX+JPHEXT,        :,:,:)
1610 ENDIF
1611 IF(LEAST_ll() .AND. .NOT. L1D) THEN
1612   XLBXUM(ILBX-NRIMX-JPHEXT+1:ILBX,:,:)     = XUT(NIU-NRIMX-JPHEXT+1:NIU,    :,:)
1613   XLBXVM(ILBX-NRIMX-JPHEXT+1:ILBX,:,:)     = XVT(NIU-NRIMX-JPHEXT+1:NIU,    :,:)
1614   XLBXWM(ILBX-NRIMX-JPHEXT+1:ILBX,:,:)     = XWT(NIU-NRIMX-JPHEXT+1:NIU,    :,:)
1615   XLBXTHM(ILBX-NRIMX-JPHEXT+1:ILBX,:,:)   = XTHT(NIU-NRIMX-JPHEXT+1:NIU,    :,:)
1616   XLBXRM(ILBX-NRIMX-JPHEXT+1:ILBX,:,:,:)   = XRT(NIU-NRIMX-JPHEXT+1:NIU,    :,:,:)
1617 ENDIF
1618 IF(LSOUTH_ll() .AND. .NOT. L1D .AND. .NOT. L2D) THEN
1619   XLBYUM(:,1:NRIMY+JPHEXT,        :)     = XUT(:,1:NRIMY+JPHEXT,      :)
1620   XLBYVM(:,1:NRIMY+JPHEXT,        :)     = XVT(:,2:NRIMY+JPHEXT+1,      :)
1621   XLBYWM(:,1:NRIMY+JPHEXT,        :)     = XWT(:,1:NRIMY+JPHEXT,  :)
1622   XLBYTHM(:,1:NRIMY+JPHEXT,        :)    = XTHT(:,1:NRIMY+JPHEXT,      :)
1623   XLBYRM(:,1:NRIMY+JPHEXT,        :,:)   = XRT(:,1:NRIMY+JPHEXT,      :,:)
1624 ENDIF
1625 IF(LNORTH_ll().AND. .NOT. L1D .AND. .NOT. L2D) THEN
1626   XLBYUM(:,ILBY-NRIMY-JPHEXT+1:ILBY,:)     = XUT(:,NJU-NRIMY-JPHEXT+1:NJU,  :)
1627   XLBYVM(:,ILBY-NRIMY-JPHEXT+1:ILBY,:)     = XVT(:,NJU-NRIMY-JPHEXT+1:NJU,  :)
1628   XLBYWM(:,ILBY-NRIMY-JPHEXT+1:ILBY,:)     = XWT(:,NJU-NRIMY-JPHEXT+1:NJU,  :)
1629   XLBYTHM(:,ILBY-NRIMY-JPHEXT+1:ILBY,:)    = XTHT(:,NJU-NRIMY-JPHEXT+1:NJU,  :)
1630   XLBYRM(:,ILBY-NRIMY-JPHEXT+1:ILBY,:,:)   = XRT(:,NJU-NRIMY-JPHEXT+1:NJU,  :,:)
1631 ENDIF
1632 DO JSV = 1, NSV
1633   IF(LWEST_ll() .AND. .NOT. L1D) &
1634   XLBXSVM(1:NRIMX+JPHEXT,        :,:,JSV)   = XSVT(1:NRIMX+JPHEXT,        :,:,JSV)
1635   IF(LEAST_ll() .AND. .NOT. L1D) &
1636   XLBXSVM(ILBX-NRIMX-JPHEXT+1:ILBX,:,:,JSV)   = XSVT(NIU-NRIMX-JPHEXT+1:NIU,    :,:,JSV)
1637   IF(LSOUTH_ll() .AND. .NOT. L1D .AND. .NOT. L2D) &
1638   XLBYSVM(:,1:NRIMY+JPHEXT,        :,JSV)   = XSVT(:,1:NRIMY+JPHEXT,      :,JSV)
1639   IF(LNORTH_ll() .AND. .NOT. L1D .AND. .NOT. L2D) &
1640   XLBYSVM(:,ILBY-NRIMY-JPHEXT+1:ILBY,:,JSV)   = XSVT(:,NJU-NRIMY-JPHEXT+1:NJU,  :,JSV)
1641 END DO
1642 !
1643 !
1644 !*       5.8   Add a perturbation to a basic state :
1645 !
1646 IF(LPERTURB) CALL SET_PERTURB(CEXPRE)
1647 !
1648 !
1649 !*       5.9   Anelastic correction and pressure:
1650 !
1651 CALL ICE_ADJUST_BIS(XPABST,XTHT,XRT)
1652 IF ( .NOT. L1D ) CALL PRESSURE_IN_PREP(XDXX,XDYY,XDZX,XDZY,XDZZ)
1653 CALL ICE_ADJUST_BIS(XPABST,XTHT,XRT)
1654 !
1655 !
1656 !*       5.10  Compute THETA, vapor and cloud mixing ratio
1657 !
1658 IF (CIDEAL == 'RSOU') THEN
1659     ALLOCATE(ZEXN(NIU,NJU,NKU))         
1660   ALLOCATE(ZT(NIU,NJU,NKU))  
1661   ALLOCATE(ZTHL(NIU,NJU,NKU))              
1662   ALLOCATE(ZRT(NIU,NJU,NKU))              
1663   ALLOCATE(ZCPH(NIU,NJU,NKU))        
1664   ALLOCATE(ZLVOCPEXN(NIU,NJU,NKU))        
1665   ALLOCATE(ZLSOCPEXN(NIU,NJU,NKU))  
1666   ALLOCATE(ZFRAC_ICE(NIU,NJU,NKU))
1667   ALLOCATE(ZRSATW(NIU,NJU,NKU))
1668   ALLOCATE(ZRSATI(NIU,NJU,NKU))             
1669   ZRT=XRT(:,:,:,1)+XRT(:,:,:,2)+XRT(:,:,:,4)
1670   ZEXN=(XPABST/XP00) ** (XRD/XCPD)
1671   ZT=XTHT*(XPABST/XP00)**(XRD/XCPD)
1672   ZCPH=XCPD+ XCPV * XRT(:,:,:,1)+ XCL *XRT(:,:,:,2)  + XCI * XRT(:,:,:,4)
1673   ZLVOCPEXN = (XLVTT + (XCPV-XCL) * (ZT-XTT))/(ZCPH*ZEXN)
1674   ZLSOCPEXN = (XLSTT + (XCPV-XCI) * (ZT-XTT))/(ZCPH*ZEXN)
1675   ZTHL=XTHT-ZLVOCPEXN*XRT(:,:,:,2)-ZLSOCPEXN*XRT(:,:,:,4)
1676   DEALLOCATE(ZEXN)         
1677   DEALLOCATE(ZT)       
1678   DEALLOCATE(ZCPH)        
1679   DEALLOCATE(ZLVOCPEXN)        
1680   DEALLOCATE(ZLSOCPEXN)
1681   CALL TH_R_FROM_THL_RT_3D('T',ZFRAC_ICE,XPABST,ZTHL,ZRT,XTHT,XRT(:,:,:,1), &
1682                             XRT(:,:,:,2),XRT(:,:,:,4),ZRSATW, ZRSATI)
1683   DEALLOCATE(ZTHL) 
1684   DEALLOCATE(ZRT)
1685 ! Coherence test
1686   IF ((.NOT. LUSERI) ) THEN
1687     IF (MAXVAL(XRT(:,:,:,4))/= 0) THEN
1688        WRITE(NLUOUT,FMT=*) "*********************************"             
1689        WRITE(NLUOUT,FMT=*) 'WARNING'      
1690        WRITE(NLUOUT,FMT=*) 'YOU HAVE LUSERI=FALSE '
1691        WRITE(NLUOUT,FMT=*) ' BUT WITH YOUR RADIOSOUNDING Ri/=0'
1692        WRITE(NLUOUT,FMT=*) MINVAL(XRT(:,:,:,4)),MAXVAL(XRT(:,:,:,4))
1693        WRITE(NLUOUT,FMT=*) "*********************************"       
1694     ENDIF  
1695   ENDIF
1696   IF ((.NOT. LUSERC)) THEN
1697     IF (MAXVAL(XRT(:,:,:,2))/= 0) THEN          
1698       WRITE(NLUOUT,FMT=*) "*********************************"
1699       WRITE(NLUOUT,FMT=*) 'WARNING'    
1700       WRITE(NLUOUT,FMT=*) 'YOU HAVE LUSERC=FALSE '
1701       WRITE(NLUOUT,FMT=*) 'BUT WITH YOUR RADIOSOUNDING RC/=0'
1702       WRITE(NLUOUT,FMT=*) MINVAL(XRT(:,:,:,2)),MAXVAL(XRT(:,:,:,2))      
1703       WRITE(NLUOUT,FMT=*) "*********************************"
1704     ENDIF  
1705   ENDIF
1706       ! on remet les bonnes valeurs pour NRR
1707   IF(CCLOUD=='NONE') NRR=1
1708   IF(CCLOUD=='REVE') NRR=2
1709 END IF
1710 !
1711 !-------------------------------------------------------------------------------
1712 !
1713 !*       6.    INITIALIZE SCALAR VARIABLES FOR CHEMISTRY
1714 !              -----------------------------------------
1715 !
1716 !  before calling chemistry
1717 CCONF = 'START'
1718 CSTORAGE_TYPE='TT'                  
1719 CALL CLOSE_ll(CEXPRE,IOSTAT=NRESP)  ! Close the EXPRE file 
1720 !
1721 IF ( LCH_INIT_FIELD ) CALL CH_INIT_FIELD_n(1, NLUOUT, NVERB)
1722 !
1723 !-------------------------------------------------------------------------------
1724 !
1725 !*       7.    WRITE THE FMFILE 
1726 !              ----------------
1727 !
1728 CALL SECOND_MNH2(ZTIME1)
1729 !
1730 NNPRAR = 22 + 2*(NRR+NSV)   &    ! 22 = number of grid variables + reference 
1731        + 8 + 17                  ! state variables + dimension variables
1732                                  ! 2*(8+NRR+NSV) + 1 = number of prognostic
1733                                  ! variables at time t and t-dt
1734 NTYPE=1
1735 CDESFM=ADJUSTL(ADJUSTR(CINIFILE)//'.des')
1736 !
1737 CALL FMOPEN_ll(CINIFILE,'WRITE',CLUOUT,NNPRAR,NTYPE,NVERB,NNINAR,NRESP)
1738 !
1739 CALL WRITE_DESFM_n(1,CDESFM,CLUOUT)
1740 !
1741 #ifdef MNH_NCWRIT
1742 NC_WRITE = LNETCDF 
1743 CALL WRITE_LFIFM_n(CINIFILE,'                            ')  ! There is no DAD model for PREP_IDEAL_CASE 
1744 IF ( LNETCDF ) THEN
1745    DEF_NC=.FALSE.
1746    CALL WRITE_LFIFM_n(CINIFILE,'                            ')  ! There is no DAD model for PREP_IDEAL_CASE 
1747    DEF_NC=.TRUE.
1748 END IF
1749 #else
1750 CALL WRITE_LFIFM_n(CINIFILE,'                            ')  ! There is no DAD model for PREP_IDEAL_CASE 
1751 #endif
1752 !
1753 CALL SECOND_MNH2(ZTIME2)
1754 !
1755 XT_STORE = XT_STORE + ZTIME2 - ZTIME1
1756 !
1757 !-------------------------------------------------------------------------------
1758 !
1759 !*     8.     EXTERNALIZED SURFACE
1760 !             --------------------
1761 !
1762 !
1763 IF (CSURF =='EXTE') THEN
1764   IF (LEN_TRIM(CINIFILEPGD)==0) THEN
1765     IF (LEN_TRIM(CPGD_FILE)/=0) THEN
1766       CINIFILEPGD=CPGD_FILE
1767     ELSE
1768       WRITE(NLUOUT,FMT=*) 'STOP : CINIFILEPGD needed in NAM_LUNITn'
1769       !callabortstop
1770       CALL CLOSE_ll(CLUOUT,IOSTAT=IRESP)
1771       CALL ABORT
1772       STOP
1773     ENDIF      
1774   ENDIF
1775   CALL ALLOC_SURFEX(1)
1776   CALL READ_ALL_NAMELISTS('MESONH','PRE',.FALSE.)                                     
1777   ! Switch to model 1 surface variables
1778   CALL GOTO_SURFEX(1,.TRUE.)
1779   !* definition of physiographic fields
1780   ! computed ...
1781   IF (LEN_TRIM(CPGD_FILE)==0 .OR. .NOT. LREAD_GROUND_PARAM) THEN
1782     CPGDFILE = CINIFILE
1783     CALL PGD_GRID_SURF_ATM('MESONH',CINIFILE,'MESONH',.TRUE.)
1784 !   CALL SPLIT_GRID('MESONH')
1785     CALL PGD_SURF_ATM     ('MESONH',CINIFILE,'MESONH',.TRUE.)
1786     CPGDFILE = CINIFILEPGD                                   
1787   ELSE
1788   ! ... or read from file.
1789     CPGDFILE = CPGD_FILE
1790     CALL INIT_PGD_SURF_ATM('MESONH','PGD',                         &
1791                             '                            ','      ',&
1792                             TDTCUR%TDATE%YEAR, TDTCUR%TDATE%MONTH,  &
1793                             TDTCUR%TDATE%DAY, TDTCUR%TIME           )
1794 !
1795   END IF
1796   !
1797   !* forces orography from atmospheric file
1798   IF (.NOT. LREAD_ZS) CALL MNHPUT_ZS_n
1799   !
1800   ! on ecrit un nouveau fichier PGD que s'il n'existe pas
1801   IF (LEN_TRIM(CPGD_FILE)==0 .OR. .NOT. LREAD_GROUND_PARAM) THEN
1802     !* writing of physiographic fields in the file
1803     CSTORAGE_TYPE='PG'
1804     COUTFMFILE = CINIFILEPGD
1805     CALL FMOPEN_ll(CINIFILEPGD,'WRITE',CLUOUT,NNPRAR,NTYPE,NVERB,NNINAR,NRESP)  
1806 #ifdef MNH_NCWRIT
1807     CALL FMWRIT(CINIFILEPGD,'PROGRAM     ',CLUOUT,'--',CPROGRAM,0,1,' ',NRESP)
1808     CALL FMWRIT(CINIFILEPGD,'SURF        ',CLUOUT,'--','EXTE',0,1,' ',NRESP)
1809     CALL FMWRIT(CINIFILEPGD,'L1D         ',CLUOUT,'--',L1D,0,1,' ',NRESP)
1810     CALL FMWRIT(CINIFILEPGD,'L2D         ',CLUOUT,'--',L2D,0,1,' ',NRESP)
1811     CALL FMWRIT(CINIFILEPGD,'PACK        ',CLUOUT,'--',LPACK,0,1,' ',NRESP)
1812     CALL WRITE_HGRID(1,CINIFILEPGD,' ')
1813     NC_FILE='sf1'
1814     NC_WRITE=LNETCDF
1815     CALL WRITE_PGD_SURF_ATM_n('MESONH')
1816     IF ( LNETCDF ) THEN
1817        DEF_NC=.FALSE.
1818        CALL WRITE_PGD_SURF_ATM_n('MESONH')
1819        DEF_NC=.TRUE.
1820        NC_WRITE = .FALSE.
1821     END IF
1822 #else
1823     CALL FMWRIT(CINIFILEPGD,'PROGRAM     ',CLUOUT,'--',CPROGRAM,0,1,' ',NRESP)
1824     CALL FMWRIT(CINIFILEPGD,'SURF        ',CLUOUT,'--','EXTE',0,1,' ',NRESP)
1825     CALL FMWRIT(CINIFILEPGD,'L1D         ',CLUOUT,'--',L1D,0,1,' ',NRESP)
1826     CALL FMWRIT(CINIFILEPGD,'L2D         ',CLUOUT,'--',L2D,0,1,' ',NRESP)
1827     CALL FMWRIT(CINIFILEPGD,'PACK        ',CLUOUT,'--',LPACK,0,1,' ',NRESP)
1828     CALL WRITE_HGRID(1,CINIFILEPGD,' ')
1829     CALL WRITE_PGD_SURF_ATM_n('MESONH')
1830 #endif
1831   CSTORAGE_TYPE='TT'
1832   ENDIF
1833   
1834   !
1835   !* deallocation of physiographic fields
1836   CALL DEALLOC_SURF_ATM_n
1837   !
1838   !* rereading of physiographic fields and definition of prognostic fields
1839   !* writing of all surface fields
1840   COUTFMFILE = CINIFILE
1841   CALL PREP_SURF_MNH('                            ','      ')
1842   CALL DEALLOC_SURFEX
1843 ELSE
1844   CSURF = "NONE"
1845 END IF
1846 !
1847 !-------------------------------------------------------------------------------
1848 !
1849 !*     9.     CLOSES THE FILE
1850 !             ---------------
1851 !
1852 IF (CSURF =='EXTE' .AND. (LEN_TRIM(CPGD_FILE)==0 .OR. .NOT. LREAD_GROUND_PARAM)) THEN
1853   CALL FMCLOS_ll(CINIFILEPGD,'KEEP',CLUOUT,NRESP)
1854 ENDIF
1855 CALL FMCLOS_ll(CINIFILE,'KEEP',CLUOUT,NRESP)
1856 IF( LEN_TRIM(CPGD_FILE) /= 0 ) THEN
1857   CALL FMCLOS_ll(CPGD_FILE,'KEEP',CLUOUT,NRESP)
1858 ENDIF
1859 !
1860 !
1861 !-------------------------------------------------------------------------------
1862 !
1863 !*      10.    PRINTS ON OUTPUT-LISTING
1864 !              ------------------------
1865 !
1866 IF (NVERB >= 5) THEN
1867   WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE: LCARTESIAN,CIDEAL,CZS=', &
1868                                     LCARTESIAN,CIDEAL,CZS 
1869   WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE: LUSERV=',LUSERV
1870   WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE: XLON0,XLAT0,XBETA,XRPK,XLONORI,XLATORI=', &
1871                                     XLON0,XLAT0,XBETA,XRPK,XLONORI,XLATORI
1872   WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE: XDELTAX,XDELTAY=',XDELTAX,XDELTAY
1873   WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE: NVERB=',NVERB
1874   IF(LCARTESIAN) THEN
1875     WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE: No map projection used.'
1876   ELSE
1877     IF (XRPK == 1.) THEN
1878       WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE: Polar stereo used.'
1879     ELSE IF (XRPK == 0.) THEN
1880       WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE: Mercator used.'
1881     ELSE
1882       WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE: Lambert used, cone factor=',XRPK 
1883     END IF
1884   END IF
1885 END IF
1886 !
1887 IF (NVERB >= 5) THEN
1888   WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE: IIB, IJB, IKB=',NIB,NJB,NKB
1889   WRITE(NLUOUT,FMT=*) 'PREP_IDEAL_CASE: IIU, IJU, IKU=',NIU,NJU,NKU
1890 END IF
1891 !
1892 !
1893 !*       28.1   print statistics!
1894 !
1895   !
1896   CALL SECOND_MNH2(ZTIME2)
1897   XT_START=XT_START+ZTIME2-ZEND
1898   !
1899   ! Set File Timing OUTPUT
1900   !
1901   CALL SET_ILUOUT_TIMING(NLUOUT)
1902   !
1903   ! Compute global time
1904   !
1905   CALL TIME_STAT_ll(XT_START,ZTOT)
1906   !
1907   !
1908   IMI = 1
1909   CALL TIME_HEADER_ll(IMI)
1910   !
1911   CALL TIME_STAT_ll(XT_STORE,ZTOT,      ' STORE-FIELDS','=')
1912   CALL  TIMING_SEPARATOR('+')
1913   CALL  TIMING_SEPARATOR('+')  
1914   WRITE(YMI,FMT="(I0)") IMI
1915   CALL TIME_STAT_ll(XT_START,ZTOT,      ' MODEL'//YMI,'+')
1916   CALL  TIMING_SEPARATOR('+')
1917   CALL  TIMING_SEPARATOR('+')
1918   CALL  TIMING_SEPARATOR('+')
1919 WRITE(NLUOUT,FMT=*) ' '
1920 WRITE(NLUOUT,FMT=*) '****************************************************'
1921 WRITE(NLUOUT,FMT=*) '* PREP_IDEAL_CASE: PREP_IDEAL_CASE ENDS CORRECTLY. *'
1922 WRITE(NLUOUT,FMT=*) '****************************************************'
1923 !
1924 CALL CLOSE_ll(CLUOUT,IOSTAT=NRESP)
1925 CALL END_PARA_ll(IINFO_ll)
1926 !
1927
1928    !callabortstop
1929    !JUAN CALL ABORT
1930 STOP
1931 !
1932 END PROGRAM PREP_IDEAL_CASE