Juan 14/04/2016: removed debug variable not allocated
[MNH-git_open_source-lfs.git] / src / MNH / spawn_zs.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 MODULE MODI_SPAWN_ZS
12 !###################
13 !
14 INTERFACE
15 !
16      SUBROUTINE SPAWN_ZS (KXOR,KXEND,KYOR,KYEND,KDXRATIO,KDYRATIO,KDIMX_C,KDIMY_C,HLBCX,HLBCY,&
17                           HLUOUT,PZS1_F,PZS2_C,HFIELD,PZS2_LS   )
18 !
19 INTEGER,   INTENT(IN)  :: KXOR,KXEND !  horizontal position (i,j) of the ORigin and END
20 INTEGER,   INTENT(IN)  :: KYOR,KYEND ! of the model 2 domain, relative to model 1
21 INTEGER,   INTENT(IN)  :: KDXRATIO   !  x and y-direction Resolution ratio
22 INTEGER,   INTENT(IN)  :: KDYRATIO   ! between model 2 and model 1
23 INTEGER,   INTENT(IN)  :: KDIMX_C    ! dimension (X dir) of local son subdomain in father grid
24 INTEGER,   INTENT(IN)  :: KDIMY_C    ! dimension (Y dir) of local son subdomain in father grid
25 CHARACTER(LEN=4),DIMENSION(2),INTENT(IN):: HLBCX, HLBCY  ! X- and Y-direc LBC
26 CHARACTER(LEN=*),     INTENT(IN)  :: HLUOUT  ! output-listing file
27 REAL, DIMENSION(:,:), INTENT(IN)  :: PZS1_F    ! model 1 orography
28 REAL, DIMENSION(:,:), INTENT(OUT) :: PZS2_C    ! interpolated orography with iterative correction
29 CHARACTER(LEN=6),     INTENT(IN)  :: HFIELD ! name of the field to nest
30 REAL, DIMENSION(:,:), INTENT(OUT),OPTIONAL  :: PZS2_LS ! interpolated orography
31 !
32 END SUBROUTINE SPAWN_ZS
33 !
34 END INTERFACE
35 !
36 END MODULE MODI_SPAWN_ZS
37 !
38 !
39 !     #########################################################################
40      SUBROUTINE SPAWN_ZS (KXOR,KXEND,KYOR,KYEND,KDXRATIO,KDYRATIO,KDIMX_C,KDIMY_C,HLBCX,HLBCY,&
41                           HLUOUT,PZS1_F,PZS2_C,HFIELD,PZS2_LS   )
42 !     #########################################################################
43 !
44 !!****  *SPAWN_ZS * - subroutine to spawn zs field
45 !!
46 !!    PURPOSE
47 !!    -------
48 !!
49 !!      This routine defines the information necessary to generate the model 2
50 !!    grid, consistently with the spawning model 1.
51 !!      The longitude and latitude of the model 2 origine are computed from
52 !!    the model 1. Then the grid in the conformal projection and terrain
53 !!    following coordinates (XHAT,YHAT and ZHAT) and orography, are interpolated
54 !!    from the model 1 grid and orography knowledge.
55 !!
56 !!**  METHOD
57 !!    ------
58 !!
59 !!      The model 2 variables are transmitted by argument (P or K prefixes),
60 !!    while the ones of model 1 are declared through calls to MODD_...
61 !!    (X or N prefixes)
62 !!
63 !!      For the case where the resolution ratio between models is 1,
64 !!    the horizontal interpolation becomes a simple equality.
65 !!      For the general case where resolution ratio is not egal to one,
66 !!    grid and orography are interpolated as follows:
67 !!         - linear interpolation for XHAT and YHAT
68 !!         - identity for ZHAT (no vertical spawning)
69 !!            2 types of interpolations can be used:
70 !!                 1. Clark and Farley (JAS 1984) on 9 points
71 !!                 2. Bikhardt on 16 points
72 !!
73 !!    EXTERNAL
74 !!    --------
75 !!
76 !!      FMLOOK        : to recover a logical unit number
77 !!      Routine BIKHARDT2     : to perform horizontal interpolations
78 !!
79 !!
80 !!    IMPLICIT ARGUMENTS
81 !!    ------------------
82 !!
83 !!      Module MODD_PARAMETERS : contains parameters
84 !!      Module MODD_CONF       : contains models configuration
85 !!      Module MODD_LUNIT2     : contains unit numbers of model 2 files
86 !!
87 !!    REFERENCE
88 !!    ---------
89 !!
90 !!       Book1 of the documentation
91 !!       PROGRAM SPAWN_ZS (Book2 of the documentation)
92 !!
93 !!    AUTHOR
94 !!    ------
95 !!
96 !!       V. Masson    * METEO-FRANCE *
97 !!
98 !!    MODIFICATIONS
99 !!    -------------
100 !!
101 !!      Original     12/01/05
102 !      Modification    20/05/06 Remove Clark and Farley interpolation
103 !      Modification    2014 M.Faivre : parallelizattion attempt
104 !      Modification    10/02/15 M. Moge : paralellization
105 !!      J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
106 !-------------------------------------------------------------------------------
107 !
108 !*       0.     DECLARATIONS
109 !               ------------
110 !
111 USE MODD_PARAMETERS, ONLY : JPHEXT       ! Declarative modules
112 USE MODD_CONF,       ONLY : NVERB
113 !
114 USE MODD_BIKHARDT_n
115 !
116 USE MODI_BIKHARDT
117 USE MODI_ZS_BOUNDARY
118 !
119 USE MODE_MODELN_HANDLER
120 !
121 USE MODE_MPPDB
122 USE MODD_VAR_ll
123 USE MODE_ll
124 USE MODD_LBC_n
125 USE MODD_NESTING
126 USE MODE_EXCHANGE_ll
127 USE MODE_EXTRAPOL
128 !
129 IMPLICIT NONE
130 !$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
131 !$20140624 renaming for VARS :
132 !  frame=Father -> _F when DIMS = IDIMX,Y
133 !  projection from |grid1 to |grid2 :
134 !  obtained with SET_LSFIELD_1WAYn + LS_FORCING 
135 !  frame=Son    -> _C when DIMS = IOR,END
136 !  projection from |grid2 to |grid1 :
137 !  obtained with SET_LSFIELD_2WAYn + LS_FEEDBACK
138 !$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
139 !
140 !*       0.1   Declarations of dummy arguments :
141 !
142 INTEGER,   INTENT(IN)  :: KXOR,KXEND !  horizontal position (i,j) of the ORigin and END
143 INTEGER,   INTENT(IN)  :: KYOR,KYEND ! of the model 2 domain, relative to model 1
144 INTEGER,   INTENT(IN)  :: KDXRATIO   !  x and y-direction Resolution ratio
145 INTEGER,   INTENT(IN)  :: KDYRATIO   ! between model 2 and model 1
146 INTEGER,   INTENT(IN)  :: KDIMX_C    ! dimension (X dir) of local son subdomain in father grid
147 INTEGER,   INTENT(IN)  :: KDIMY_C    ! dimension (Y dir) of local son subdomain in father grid
148 CHARACTER(LEN=4),DIMENSION(2),INTENT(IN):: HLBCX, HLBCY  ! X- and Y-direc LBC
149 CHARACTER(LEN=*),     INTENT(IN)  :: HLUOUT  ! output-listing file
150 REAL, DIMENSION(:,:), INTENT(IN)  :: PZS1_F    ! model 1 orography
151 REAL, DIMENSION(:,:), INTENT(OUT) :: PZS2_C    ! interpolated orography with iterative correction
152 CHARACTER(LEN=6),     INTENT(IN)  :: HFIELD ! name of the field to nest
153 REAL, DIMENSION(:,:), INTENT(OUT),OPTIONAL  :: PZS2_LS ! interpolated orography
154 !
155 !*       0.2    Declarations of local variables for print on FM file
156 !
157 INTEGER :: ILUOUT   ! Logical unit number for the output listing
158 INTEGER :: IRESP    ! Return codes in FM routines
159 !
160 REAL, DIMENSION(:,:), ALLOCATABLE :: ZZS2_LS ! interpolated orography
161 REAL, DIMENSION(:,:), ALLOCATABLE :: PZS1_C ! model 1 orography resticted to the grid of model 2
162 REAL, DIMENSION(:,:), ALLOCATABLE :: ZZS1_C ! zs of model 1 at iteration n or n+1 in GRID2
163 REAL, DIMENSION(:,:), ALLOCATABLE :: ZZS2_C ! averaged zs of model 2 at iteration n
164 REAL, DIMENSION(:,:), ALLOCATABLE :: ZDZS_C ! difference between PZS1 and ZZS2
165 !$20140617 ZTZS1 result of SET_LSFIELD_1WAY_ll(PZS1)
166 !JUAN REAL, DIMENSION(:,:), ALLOCATABLE :: ZTZS1_C
167 !$20140704 ZDZS_3D to use MAX_ll(array3D arg)
168 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZDZS_3D
169 !$
170 INTEGER                :: IXMIN, IXMAX ! indices to interpolate the
171                                     ! modified orography on model 2
172                                     ! domain to model 1 grid
173 !JUAN A REVOIR TODO_JPHEXT /!\ /!\
174 ! <<<<<<< spawn_zs.f90
175 INTEGER                :: JI,JEPSX  ! Loop index in x direction
176 INTEGER                :: JJ,JEPSY  ! Loop index in y direction
177 INTEGER                :: JCOUNTER  ! counter for iterative method
178 REAL                   :: ZRELAX             ! relaxation factor
179 INTEGER                :: JMAXITER = 2000    ! maximum number of iterations
180 !
181 INTEGER, DIMENSION(2)  :: IZSMAX
182 INTEGER                :: IMI     ! current model index
183 !$20140604
184 INTEGER                :: KMI,IDIMX_C,IDIMY_C
185 !$20140602
186 INTEGER                :: PZS1_FSIZE1
187 INTEGER                :: PZS1_FSIZE2
188 !$20140603
189 INTEGER                :: IINFO_ll
190 !$20140619
191 TYPE(LIST_ll), POINTER :: TZFIELDS_ll => NULL()  ! list of fields to exchange
192 !$20140623
193 INTEGER                :: IXOR_F,IXEND_F
194 INTEGER                :: IYOR_F,IYEND_F
195 INTEGER                :: KDXRATIO_C, KDYRATIO_C
196 !$20140704
197 !$20140711 not INT, REAL !!
198 REAL                   :: ZMAXVAL
199 REAL                   :: LOCMAXVAL
200 !$20140801
201 INTEGER                :: IORX, IORY, IIBINT,IJBINT,IIEINT,IJEINT
202 INTEGER                :: IXOR_C_ll, IXEND_C_ll  ! origin and end of the local subdomain of the child model 2
203 INTEGER                :: IYOR_C_ll, IYEND_C_ll  ! relative to the father model 1
204 INTEGER                :: IINFO  ! return code of // routines
205 !
206 TYPE(LIST_ll), POINTER :: TZZSFIELD_ll   ! list of fields to exchange
207 TYPE(HALO2LIST_ll), POINTER :: TZZSHALO2_ll   ! needed for update_halo2_ll
208
209 REAL, DIMENSION(:,:), ALLOCATABLE :: ZZS1CHILDGRID_C  ! copy of ZZS1_C extended to the whole child domain
210 INTEGER                :: JI2INF,JI2SUP
211 INTEGER                :: JJ2INF,JJ2SUP
212 INTEGER               :: IXSIZE,IYSIZE
213 !
214 KDXRATIO_C=KDXRATIO
215 KDYRATIO_C=KDYRATIO
216 ZMAXVAL=1000.
217 !$
218 !!$=======
219 !!$INTEGER             :: JI,JEPSX  ! Loop index in x direction
220 !!$INTEGER             :: JJ,JEPSY  ! Loop index in y direction
221 !!$INTEGER             :: JCOUNTER  ! counter for iterative method
222 !!$REAL                :: ZRELAX             ! relaxation factor
223 !!$INTEGER             :: JMAXITER = 2000    ! maximum number of iterations
224 !!$!
225 !!$INTEGER, DIMENSION(2) :: IZSMAX
226 !!$INTEGER               :: IMI     ! current model index
227 !!$!
228 !!$INTEGER               :: IXSIZE,IYSIZE
229 !!$INTEGER                          :: INFO_ll                ! error return code
230 !!$>>>>>>> 1.1.4.1.18.2.2.1
231 !-------------------------------------------------------------------------------
232 !
233 !*       1.    PROLOGUE:
234 !              ---------
235 !
236 IMI = GET_CURRENT_MODEL_INDEX()
237 CALL GOTO_MODEL(IMI)
238 !
239 !
240 !*       1.2  recovers logical unit number of output listing
241 !
242 CALL FMLOOK_ll(HLUOUT,HLUOUT,ILUOUT,IRESP)
243 !
244 !-------------------------------------------------------------------------------
245 !
246 !*       2.    COMPUTATION:
247 !              ---------
248 !
249 !*       2.1   Purely interpolated zs:
250 !              -----------------------
251 !
252 ALLOCATE(ZZS2_LS(SIZE(PZS2_C,1),SIZE(PZS2_C,2)))
253 PZS1_FSIZE1=SIZE(PZS1_F,1)
254 PZS1_FSIZE2=SIZE(PZS1_F,2)
255 !
256 !
257 !
258 ! This is one way to do it, but then the compute load are not well balanced
259 ! Each process computes the interpolation on the intersection of the global
260 ! child model with its part of the father model
261 !CALL BIKHARDT (XBMX1,XBMX2,XBMX3,XBMX4,XBMY1,XBMY2,XBMY3,XBMY4, &
262 !               XBFX1,XBFX2,XBFX3,XBFX4,XBFY1,XBFY2,XBFY3,XBFY4, &
263 !               KXOR,KYOR,KXEND,KYEND,KDXRATIO,KDYRATIO,1,       &
264 !               HLBCX,HLBCY,PZS1_F,ZZS2_LS)
265 !
266 ! We want instead that each process computes the interpolation for its part
267 ! of the child model
268 ! Then we have to communicate the values of PZS1_F on each subdomain of the
269 ! child model to the corresponding process
270 ! before calling BIKHARDT on the local subdomain of the child model
271 !
272 !*      1 GATHER LS FIELD FOR THE CHILD MODEL KMI
273 !       1.1  Must be on the father model to call get_child_dim
274 CALL GOTO_MODEL(NDAD(IMI))
275 !$20140623 KMI is DAD, IMI=son !!
276 !$20140623 use IMI not KMI
277 CALL GO_TOMODEL_ll(NDAD(IMI), IINFO_ll)
278 IDIMX_C = KDIMX_C! + 2*(JPHEXT+1) !KXEND-KXOR-1
279 IDIMY_C = KDIMY_C! + 2*(JPHEXT+1) !KYEND-KYOR-1
280 !CALL GET_CHILD_DIM_ll(IMI, IDIMX_C, IDIMY_C, IINFO_ll)
281 !
282 !         1.3  Specify the ls "source" fields and receiver fields
283 !
284 ALLOCATE(ZZS1_C(IDIMX_C,IDIMY_C))
285 ZZS1_C(:,:)=0.
286 CALL SET_LSFIELD_1WAY_ll(PZS1_F, ZZS1_C, IMI)
287 CALL MPPDB_CHECK2D(PZS1_F,"SPAWN_ZS:PZS1_F",PRECISION)
288 !        1.4  Communication
289 CALL LS_FORCING_ll(IMI, IINFO_ll, .TRUE.)
290 !        1.5  Back to the (current) child model
291 CALL GO_TOMODEL_ll(IMI, IINFO_ll)
292 CALL GOTO_MODEL(IMI)
293 CALL UNSET_LSFIELD_1WAY_ll()
294 !
295 !if the child grid is the whole father grid, we first need to extrapolate
296 !the data on a "pseudo halo" before doing BIKHARDT interpolation
297 !CALL EXTRAPOL_ON_PSEUDO_HALO(ZZS1_C)
298 ! <<<<<<< spawn_zs.f90
299 CALL BIKHARDT (XBMX1,XBMX2,XBMX3,XBMX4,XBMY1,XBMY2,XBMY3,XBMY4, &
300                XBFX1,XBFX2,XBFX3,XBFX4,XBFY1,XBFY2,XBFY3,XBFY4, &
301                2,2,IDIMX_C-1,IDIMY_C-1,KDXRATIO_C,KDYRATIO_C,1, &
302                HLBCX,HLBCY,ZZS1_C,ZZS2_LS)
303 !!$=======
304 !
305 !!$  CALL BIKHARDT (XBMX1,XBMX2,XBMX3,XBMX4,XBMY1,XBMY2,XBMY3,XBMY4, &
306 !!$                 XBFX1,XBFX2,XBFX3,XBFX4,XBFY1,XBFY2,XBFY3,XBFY4, &
307 !!$                 KXOR,KYOR,KXEND,KYEND,KDXRATIO,KDYRATIO,1,       &
308 !!$                 HLBCX,HLBCY,PZS1,ZZS2_LS)
309 !!$>>>>>>> 1.1.4.1.18.2.2.1
310   CALL MPPDB_CHECK2D(ZZS2_LS,"SPAWN_ZS::ZZS2_LS",PRECISION)
311 !
312 !*       4.2   New zs:
313 !              -------
314 !
315 !* we use an iterative method to insure the equality between large-scale
316 !  orography and the average of fine orography, and to be sure not to have
317 !  spurious cliffs near the coast (multiplication by xland during the
318 !  iterative process).
319 !
320 IF (KDXRATIO/=1 .OR. KDYRATIO/=1) THEN
321 !
322 !* allocations
323 ! <<<<<<< spawn_zs.f90
324    IXSIZE = IDIMX_C-2*(JPHEXT+1)
325    IYSIZE = IDIMY_C-2*(JPHEXT+1)
326    ALLOCATE(ZZS2_C(IXSIZE,IYSIZE))
327    ALLOCATE(ZDZS_C(IXSIZE,IYSIZE))
328    ALLOCATE(PZS1_C(IXSIZE,IYSIZE))
329 !!$=======
330 !!$  IXSIZE = KXEND-KXOR - 2*JPHEXT + 1
331 !!$  IYSIZE = KYEND-KYOR - 2*JPHEXT + 1
332 !!$>>>>>>> 1.1.4.1.18.2.2.1
333 !
334 !* constants
335 !
336   ZRELAX=16./13.     ! best relaxation for infinite aspect ratio.
337                      ! for dx=2, one should take 32./27. !!!
338 !
339 !* initializations of initial state
340 !
341   JCOUNTER=0
342   PZS1_C(:,:) = ZZS1_C(JPHEXT+2:IDIMX_C-JPHEXT-1,JPHEXT+2:IDIMY_C-JPHEXT-1)
343   PZS2_C=0.
344   CALL MPPDB_CHECK2D(PZS2_C,"SPAWN_ZSbefBKAT:PZS2",PRECISION)
345 !
346 !* iterative loop
347 !
348   DO
349 !
350 !* interpolation
351 !
352 ! <<<<<<< spawn_zs.f90
353     CALL BIKHARDT (XBMX1,XBMX2,XBMX3,XBMX4,XBMY1,XBMY2,XBMY3,XBMY4,  &
354                    XBFX1,XBFX2,XBFX3,XBFX4,XBFY1,XBFY2,XBFY3,XBFY4,  &
355                    2,2,IDIMX_C-1,IDIMY_C-1,KDXRATIO_C,KDYRATIO_C,1,  &
356                    HLBCX,HLBCY,ZZS1_C,PZS2_C)
357 !!$=======
358 !!$      CALL BIKHARDT (XBMX1,XBMX2,XBMX3,XBMX4,XBMY1,XBMY2,XBMY3,XBMY4, &
359 !!$                     XBFX1,XBFX2,XBFX3,XBFX4,XBFY1,XBFY2,XBFY3,XBFY4, &
360 !!$                     KXOR,KYOR,KXEND,KYEND,KDXRATIO,KDYRATIO,1,       &
361 !!$                     HLBCX,HLBCY,ZZS1,PZS2)
362       CALL MPPDB_CHECK2D(PZS2_C,"SPAWN_ZS::PZS2_C/LOOP",PRECISION)
363 !!$>>>>>>> 1.1.4.1.18.2.2.1
364     JCOUNTER=JCOUNTER+1
365 !
366 !* if orography is positive, it stays positive
367 !
368 ! <<<<<<< spawn_zs.f90
369      DO JI=1,IXSIZE
370        DO JJ=1,IYSIZE
371          IF (PZS1_C(JI,JJ)>-1.E-15) THEN
372             JI2INF = (JI-1)*KDXRATIO_C+1+JPHEXT
373             JI2SUP = JI*KDXRATIO_C+JPHEXT
374             JJ2INF = (JJ-1)*KDYRATIO_C+1+JPHEXT
375             JJ2SUP = JJ*KDYRATIO_C+JPHEXT
376             PZS2_C(JI2INF:JI2SUP,JJ2INF:JJ2SUP)= MAX(PZS2_C(JI2INF:JI2SUP,JJ2INF:JJ2SUP),0.)
377          ENDIF
378        END DO
379      END DO
380 !!$=======
381 !!$    DO JI=1,IXSIZE
382 !!$      DO JJ=1,IYSIZE
383 !!$        IF (PZS1(JI+KXOR,JJ+KYOR)>-1.E-15)                            &
384 !!$          PZS2((JI-1)*KDXRATIO+1+JPHEXT:JI*KDXRATIO+JPHEXT,             &
385 !!$               (JJ-1)*KDYRATIO+1+JPHEXT:JJ*KDYRATIO+JPHEXT)  =          &
386 !!$            MAX( PZS2((JI-1)*KDXRATIO+1+JPHEXT:JI*KDXRATIO+JPHEXT,      &
387 !!$                      (JJ-1)*KDYRATIO+1+JPHEXT:JJ*KDYRATIO+JPHEXT), 0.)
388 !!$      END DO
389 !!$    END DO
390     CALL MPPDB_CHECK2D(PZS2_C,"SPAWN_ZS::PZS2_C/POS",PRECISION)
391 !!$>>>>>>> 1.1.4.1.18.2.2.1
392 !
393 !* computation of new averaged orography
394 ! <<<<<<< spawn_zs.f90
395    ZZS2_C(:,:) = 0.
396    DO JI=1,IXSIZE
397      DO JJ=1,IYSIZE
398        DO JEPSX = (JI-1)*KDXRATIO_C+1+JPHEXT, JI*KDXRATIO_C+JPHEXT
399          DO JEPSY = (JJ-1)*KDYRATIO_C+1+JPHEXT, JJ*KDYRATIO_C+JPHEXT
400            ZZS2_C(JI,JJ) = ZZS2_C(JI,JJ) + PZS2_C(JEPSX,JEPSY)
401          END DO
402        END DO
403      END DO
404 !!$=======
405 !!$>>>>>>> 1.1.4.1.18.2.2.1
406     END DO
407 ! <<<<<<< spawn_zs.f90
408     ZZS2_C(:,:) = ZZS2_C(:,:) / (KDXRATIO_C*KDYRATIO_C)
409     !
410     ZDZS_C(:,:)=PZS1_C(:,:)-ZZS2_C(:,:)
411 !!$=======
412 !!$    ZZS2(:,:) = ZZS2(:,:) / (KDXRATIO*KDYRATIO)
413 !!$    !
414 !!$    ZDZS(:,:)=PZS1(KXOR+JPHEXT:KXEND-JPHEXT,KYOR+JPHEXT:KYEND-JPHEXT)-ZZS2(:,:)
415 !!$>>>>>>> 1.1.4.1.18.2.2.1
416 !
417 !* test to end the iterative process
418 !
419     ALLOCATE(ZDZS_3D(SIZE(ZDZS_C,1),SIZE(ZDZS_C,2),1))  ! WARNING : this is highly inefficient, this copy is unecessary
420     ZDZS_3D(:,:,1)=ZDZS_C(:,:)                          ! We could write a function MAX2D_ll or use a POINTER for ZDZS_3D
421     LOCMAXVAL=MAXVAL(ABS(ZDZS_C))
422     CALL MPI_ALLREDUCE(LOCMAXVAL,ZMAXVAL,1,MPI_PRECISION,MPI_MAX,NMNH_COMM_WORLD,IINFO_ll)
423     IF (ZMAXVAL<1.E-3) THEN
424       EXIT
425     ENDIF
426 !
427     IF (JCOUNTER>=JMAXITER) THEN
428       WRITE(ILUOUT,FMT=*) 'SPAWN_ZS: convergence of ',TRIM(HFIELD), &
429                           ' NOT obtained after',JCOUNTER,' iterations'
430       WRITE(ILUOUT,FMT=*) TRIM(HFIELD),                             &
431          ' is modified to insure egality of large scale and averaged fine field'
432 ! <<<<<<< spawn_zs.f90
433       DO JI=1,IXSIZE
434         DO JJ=1,IYSIZE
435           DO JEPSX = (JI-1)*KDXRATIO_C+1+JPHEXT, JI*KDXRATIO_C+JPHEXT
436             DO JEPSY = (JJ-1)*KDYRATIO_C+1+JPHEXT, JJ*KDYRATIO_C+JPHEXT
437               PZS2_C(JEPSX,JEPSY) = PZS2_C(JEPSX,JEPSY) + ZDZS_C(JI,JJ)
438 !!$>>>>>>> 1.1.4.1.18.2.2.1
439             END DO
440           END DO
441         END DO
442       END DO
443       !
444       EXIT
445     END IF
446 !
447 !* prints
448 !
449     IF (NVERB >=7) THEN
450       IZSMAX=MAXLOC(ABS(ZDZS_C(:,:)))
451       IF (MOD(JCOUNTER,500)==1) THEN
452         WRITE(ILUOUT,FMT='(A4,1X,A4,1X,A2,1X,A2,1X,A12,1X,A12,1X,A12)')      &
453                           'n IT','nDIV','I1','J1','   ZS1','   ZS2','   DZS'
454         WRITE(ILUOUT,FMT='(I4,1X,I4,1X,I2,1X,I2,1X,F12.7,1X,F12.7,1X,F12.7)')&
455                             JCOUNTER,COUNT(ABS(ZDZS_C(:,:))>=1.E-3),          &
456                             IZSMAX(1)+2,IZSMAX(2)+2,                          &
457 !                            PZS1(KXOR+IZSMAX(1),KYOR+IZSMAX(2)),              &
458 !                            ZTZS1_C(2+IZSMAX(1),2+IZSMAX(2)),                 &
459                             ZZS2_C(IZSMAX(1),IZSMAX(2)),                      &
460                             ZDZS_C(IZSMAX(1),IZSMAX(2))
461       ENDIF
462     END IF
463 !
464 !* correction of coarse orography
465 !
466 ! <<<<<<< spawn_zs.f90
467     ZZS1_C(JPHEXT+2:IDIMX_C-JPHEXT-1,JPHEXT+2:IDIMY_C-JPHEXT-1) = &
468     ZZS1_C(JPHEXT+2:IDIMX_C-JPHEXT-1,JPHEXT+2:IDIMY_C-JPHEXT-1) + ZRELAX * ZDZS_C(:,:)
469
470     ! update the Halo
471     ! UPDATE_HALO_ll routines only work with fields of the size of the subdomain
472     ! so we have to copy the values we want to update in a temporary field ZZS1CHILDGRID_C
473     ALLOCATE(ZZS1CHILDGRID_C(SIZE(PZS2_C,1)+2,SIZE(PZS2_C,2)+2))
474     ! TODO : renommer ZZS1CHILDGRID_C avec un nom plus explicite
475     ZZS1CHILDGRID_C = 0.
476     ! west boundary of ZZS1_C
477     DO JI=1,JPHEXT+1
478       DO JJ=1,IDIMY_C
479         ZZS1CHILDGRID_C(JI,JJ) = ZZS1_C(JI,JJ)  ! distant value, not on local physical domain
480         ZZS1CHILDGRID_C(JI+JPHEXT+1,JJ) = ZZS1_C(JI+JPHEXT+1,JJ) ! local value, on local physical domain
481       END DO
482     END DO
483     ! east boundary of ZZS1_C
484     DO JI=1,JPHEXT+1
485       DO JJ=1,IDIMY_C
486         ZZS1CHILDGRID_C(SIZE(PZS2_C,1)+2-JI+1,JJ) = ZZS1_C(IDIMX_C-JI+1,JJ)  ! distant value, not on local physical domain
487         ZZS1CHILDGRID_C(SIZE(PZS2_C,1)+2-JI+1-JPHEXT-1,JJ) = ZZS1_C(IDIMX_C-JI+1-JPHEXT-1,JJ) ! local value, on local physical domain
488       END DO
489     END DO
490     ! south boundary of ZZS1_C
491     DO JI=1,IDIMX_C
492       DO JJ=1,JPHEXT+1
493         ZZS1CHILDGRID_C(JI,JJ) = ZZS1_C(JI,JJ)  ! distant value, not on local physical domain
494         ZZS1CHILDGRID_C(JI,JJ+JPHEXT+1) = ZZS1_C(JI,JJ+JPHEXT+1) ! local value, on local physical domain
495       END DO
496     END DO
497     ! north boundary of ZZS1_C
498     DO JI=1,IDIMX_C
499       DO JJ=1,JPHEXT+1
500         ZZS1CHILDGRID_C(JI,SIZE(PZS2_C,2)+2-JJ+1) = ZZS1_C(JI,IDIMY_C-JJ+1)  ! distant value, not on local physical domain
501         ZZS1CHILDGRID_C(JI,SIZE(PZS2_C,2)+2-JJ+1-JPHEXT-1) = ZZS1_C(JI,IDIMY_C-JJ+1-JPHEXT-1) ! local value, on local physical domain
502       END DO
503     END DO
504     ! If we leave the north-east corner with zero values, UPDATE_HALO_EXTENDED_ll will
505     ! cause errors on the south-east and north-west internal border of the neigbouring processes
506     DO JI=1,JPHEXT+1
507       DO JJ=1,JPHEXT+1
508         ZZS1CHILDGRID_C(SIZE(PZS2_C,1)+2-JI+1-JPHEXT-1,SIZE(PZS2_C,2)+2-JJ+1-JPHEXT-1) &
509          = ZZS1_C(IDIMX_C-JI+1-JPHEXT-1,IDIMY_C-JJ+1-JPHEXT-1) ! local value, on local physical domain
510       END DO
511     END DO
512     !
513     NULLIFY(TZZSFIELD_ll)
514     CALL ADD2DFIELD_ll(TZZSFIELD_ll, ZZS1CHILDGRID_C)
515     CALL UPDATE_HALO_EXTENDED_ll(TZZSFIELD_ll,IINFO)
516     CALL CLEANLIST_ll(TZZSFIELD_ll)
517     ! west and east boundaries - distant points
518     DO JI=1,JPHEXT+1
519       DO JJ=JPHEXT+1,IDIMY_C-JPHEXT+1
520         ZZS1_C(JI,JJ) = ZZS1CHILDGRID_C(JI,JJ)
521         ZZS1_C(IDIMX_C-JI+1,JJ) = ZZS1CHILDGRID_C(SIZE(PZS2_C,1)+2-JI+1,JJ)
522       END DO
523     END DO
524     ! north and south boundaries - distant points
525     DO JI=JPHEXT+1,IDIMX_C-JPHEXT+1
526       DO JJ=1,JPHEXT+1
527         ZZS1_C(JI,JJ) = ZZS1CHILDGRID_C(JI,JJ)
528         ZZS1_C(JI,IDIMY_C-JJ+1) = ZZS1CHILDGRID_C(JI,SIZE(PZS2_C,2)+2-JJ+1)
529       END DO
530     END DO
531     ! "corner" halo
532     DO JI=1,JPHEXT+1
533       DO JJ=1,JPHEXT+1
534         ZZS1_C(JI,JJ) = ZZS1CHILDGRID_C(JI,JJ)
535         ZZS1_C(IDIMX_C-JI+1,JJ) = ZZS1CHILDGRID_C(SIZE(PZS2_C,1)+2-JI+1,JJ)
536         ZZS1_C(JI,IDIMY_C-JJ+1) = ZZS1CHILDGRID_C(JI,SIZE(PZS2_C,2)+2-JJ+1)
537         ZZS1_C(IDIMX_C-JI+1,IDIMY_C-JJ+1) = ZZS1CHILDGRID_C(SIZE(PZS2_C,1)+2-JI+1,SIZE(PZS2_C,2)+2-JJ+1)
538       END DO
539     END DO
540     ! corner points - distant points
541     ! we have to treat the halo points in the corner separately to have correct values
542     ! in the intersection of the halos (points (1,1), (1,2), (2,1), (2,2), (IDIMX_C,IDIMY_C), etc.)
543 !!$=======
544 !!$    ZZS1(KXOR+JPHEXT:KXEND-JPHEXT,KYOR+JPHEXT:KYEND-JPHEXT) =                             &
545 !!$         ZZS1(KXOR+JPHEXT:KXEND-JPHEXT,KYOR+JPHEXT:KYEND-JPHEXT) + ZRELAX * ZDZS(:,:)
546 !!$>>>>>>> 1.1.4.1.18.2.2.1
547     !
548     ! extrapolations (X direction)
549     !
550 ! <<<<<<< spawn_zs.f90
551     ! TODO: utiliser JPHEXT dans une boucle pour generaliser au cas ou le halo est plus grand que 1
552     IF(KXOR==1 .AND. KXEND==SIZE(PZS1_F,1) .AND. HLBCX(1)=='CYCL' ) THEN
553       !c'est pris en compte et deja fait par UPDATE_HALO_ll et UPDATE_HALO2_ll ? --------> NON
554 !!$=======
555 !!$    IF(KXOR==1 .AND. KXEND==SIZE(PZS1,1) .AND. HLBCX(1)=='CYCL' ) THEN
556 !!$      ZZS1(KXOR,KYOR+JPHEXT:KYEND-JPHEXT)  = ZZS1(KXEND-JPHEXT,KYOR+JPHEXT:KYEND-JPHEXT)
557 !!$      ZZS1(KXEND,KYOR+JPHEXT:KYEND-JPHEXT) = ZZS1(KXOR+JPHEXT,KYOR+JPHEXT:KYEND-JPHEXT)
558 !!$>>>>>>> 1.1.4.1.18.2.2.1
559     ELSE
560 ! <<<<<<< spawn_zs.f90
561       IF ( LWEST_ll() ) THEN
562         ZZS1_C(1+JPHEXT,1:IDIMY_C) = 2. * ZZS1_C(2+JPHEXT,1:IDIMY_C)  - ZZS1_C(3+JPHEXT,1:IDIMY_C)
563         ZZS1_C(JPHEXT,1:IDIMY_C)   = 2. * ZZS1_C(1+JPHEXT,1:IDIMY_C)  - ZZS1_C(2+JPHEXT,1:IDIMY_C)
564       ENDIF
565       IF ( LEAST_ll() ) THEN
566           ZZS1_C(IDIMX_C-JPHEXT,1:IDIMY_C)   = 2. * ZZS1_C(IDIMX_C-JPHEXT-1,1:IDIMY_C) - ZZS1_C(IDIMX_C-JPHEXT-2,1:IDIMY_C)
567           ZZS1_C(IDIMX_C-JPHEXT+1,1:IDIMY_C) = 2. * ZZS1_C(IDIMX_C-JPHEXT,1:IDIMY_C)   - ZZS1_C(IDIMX_C-JPHEXT-1,1:IDIMY_C)
568       ENDIF
569 !!$=======
570 !!$      ZZS1(KXOR+JPHEXT-1,KYOR+JPHEXT:KYEND-JPHEXT) =                                       &
571 !!$        2. * ZZS1(KXOR+JPHEXT,KYOR+JPHEXT:KYEND-JPHEXT)  - ZZS1(KXOR+JPHEXT+1,KYOR+JPHEXT:KYEND-JPHEXT)
572 !!$      IF(KXOR>1)                                                        &
573 !!$      ZZS1(KXOR+JPHEXT-2,KYOR+JPHEXT:KYEND-JPHEXT) =                                     &
574 !!$        2. * ZZS1(KXOR+JPHEXT-1,KYOR+JPHEXT:KYEND-JPHEXT)  - ZZS1(KXOR+JPHEXT,KYOR+JPHEXT:KYEND-JPHEXT)
575 !!$      ZZS1(KXEND-JPHEXT+1,KYOR+JPHEXT:KYEND-JPHEXT) =                                      &
576 !!$        2. * ZZS1(KXEND-JPHEXT,KYOR+JPHEXT:KYEND-JPHEXT) - ZZS1(KXEND-JPHEXT-1,KYOR+JPHEXT:KYEND-JPHEXT)
577 !!$      IF(KXEND<SIZE(PZS1,1))                                             &
578 !!$      ZZS1(KXEND-JPHEXT+2,KYOR+JPHEXT:KYEND-JPHEXT) =                                    &
579 !!$        2. * ZZS1(KXEND-JPHEXT+1,KYOR+JPHEXT:KYEND-JPHEXT) - ZZS1(KXEND-JPHEXT,KYOR+JPHEXT:KYEND-JPHEXT)
580 !!$>>>>>>> 1.1.4.1.18.2.2.1
581     END IF
582     !
583     ! extrapolations (Y direction)
584     !
585 ! <<<<<<< spawn_zs.f90
586 !    IXMIN=MAX(KXOR-1,1)
587 !    IXMAX=MIN(KXEND+1,SIZE(PZS1_F,1))
588     IF(KYOR==1 .AND. KYEND==SIZE(PZS1_F,2) .AND. HLBCY(1)=='CYCL' ) THEN
589       !c'est pris en compte et deja fait par UPDATE_HALO_ll et UPDATE_HALO2_ll ? --------> NON
590 !!$=======
591 !!$    IXMIN=MAX(KXOR-1,1)
592 !!$    IXMAX=MIN(KXEND+1,SIZE(PZS1,1))
593 !!$    IF(KYOR==1 .AND. KYEND==SIZE(PZS1,2) .AND. HLBCY(1)=='CYCL' ) THEN
594 !!$      ZZS1(IXMIN:IXMAX,KYOR)  = ZZS1(IXMIN:IXMAX,KYEND-JPHEXT)
595 !!$      ZZS1(IXMIN:IXMAX,KYEND) = ZZS1(IXMIN:IXMAX,KYOR+JPHEXT)
596 !!$>>>>>>> 1.1.4.1.18.2.2.1
597     ELSE
598 ! <<<<<<< spawn_zs.f90
599       IF ( LSOUTH_ll() ) THEN
600         ZZS1_C(1:IDIMX_C,1+JPHEXT) = 2. * ZZS1_C(1:IDIMX_C,2+JPHEXT)  - ZZS1_C(1:IDIMX_C,3+JPHEXT)
601         ZZS1_C(1:IDIMX_C,JPHEXT)   = 2. * ZZS1_C(1:IDIMX_C,1+JPHEXT)  - ZZS1_C(1:IDIMX_C,2+JPHEXT)
602       ENDIF
603       IF ( LNORTH_ll() ) THEN
604         ZZS1_C(1:IDIMX_C,IDIMY_C-JPHEXT)   = 2. * ZZS1_C(1:IDIMX_C,IDIMY_C-JPHEXT-1) - ZZS1_C(1:IDIMX_C,IDIMY_C-JPHEXT-2)
605         ZZS1_C(1:IDIMX_C,IDIMY_C-JPHEXT+1) = 2. * ZZS1_C(1:IDIMX_C,IDIMY_C-JPHEXT)   - ZZS1_C(1:IDIMX_C,IDIMY_C-JPHEXT-1)
606       ENDIF
607 !!$=======
608 !!$      ZZS1(IXMIN:IXMAX,KYOR+JPHEXT-1) =                                       &
609 !!$        2. * ZZS1(IXMIN:IXMAX,KYOR+JPHEXT)  - ZZS1(IXMIN:IXMAX,KYOR+JPHEXT+1)
610 !!$      IF(KYOR>1)                                                     &
611 !!$      ZZS1(IXMIN:IXMAX,KYOR+JPHEXT-2) =                                     &
612 !!$        2. * ZZS1(IXMIN:IXMAX,KYOR+JPHEXT-1)    - ZZS1(IXMIN:IXMAX,KYOR+JPHEXT)
613 !!$      ZZS1(IXMIN:IXMAX,KYEND-JPHEXT+1) =                                      &
614 !!$        2. * ZZS1(IXMIN:IXMAX,KYEND-JPHEXT) - ZZS1(IXMIN:IXMAX,KYEND-JPHEXT-1)
615 !!$      IF(KYEND<SIZE(PZS1,2))                                          &
616 !!$      ZZS1(IXMIN:IXMAX,KYEND-JPHEXT+2) =                                    &
617 !!$        2. * ZZS1(IXMIN:IXMAX,KYEND-JPHEXT+1)   - ZZS1(IXMIN:IXMAX,KYEND-JPHEXT)
618 !!$>>>>>>> 1.1.4.1.18.2.2.1
619     END IF
620 !
621     DEALLOCATE(ZZS1CHILDGRID_C)
622     DEALLOCATE(ZDZS_3D)
623 !
624   END DO
625 !
626   CALL ZS_BOUNDARY(PZS2_C,ZZS2_LS)
627   JI=0
628   CALL MPPDB_CHECK2D(PZS2_C,"SPAWN_ZSend:PZS2",PRECISION)
629 !
630   WRITE(ILUOUT,FMT=*) 'convergence of ',TRIM(HFIELD),' obtained after ', &
631                       JCOUNTER,' iterations'
632 !
633   DEALLOCATE(ZZS2_C)
634   DEALLOCATE(ZDZS_C)
635   DEALLOCATE(ZZS1_C)
636 END IF
637 !
638 IF (PRESENT(PZS2_LS)) PZS2_LS(:,:)=ZZS2_LS(:,:)
639 DEALLOCATE(ZZS2_LS)
640 !
641 CALL GOTO_MODEL(IMI)
642 CALL GO_TOMODEL_ll(IMI,IINFO_ll)
643 !-------------------------------------------------------------------------------
644 END SUBROUTINE SPAWN_ZS
645 !