New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
crsdom.F90 in NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS – NEMO

source: NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdom.F90 @ 10234

Last change on this file since 10234 was 10234, checked in by cbricaud, 4 years ago

correc a bug for MPP decompostion for coarsened grid for limited area model

  • Property svn:keywords set to Id
File size: 77.5 KB
Line 
1MODULE crsdom
2   !!===================================================================
3   !!                  ***  crs.F90 ***
4   !!  Purpose: Interface for calculating quantities from a 
5   !!           higher-resolution grid for the coarse grid.
6   !!
7   !!  Method:
8
9   !!  References:  Aumont, O., J.C. Orr, D. Jamous, P. Monfray
10   !!               O. Marti and G. Madec, 1998. A degradation
11   !!               approach to accelerate simulations to steady-state
12   !!               in a 3-D tracer transport model of the global ocean.
13   !!               Climate Dynamics, 14:101-116.
14   !!  History:
15   !!       Original.   May 2012.  (J. Simeon, C. Calone, G. Madec, C. Ethe)
16   !!===================================================================
17
18   USE dom_oce        ! ocean space and time domain and to get jperio
19   USE wrk_nemo       ! work arrays
20   USE crs            ! domain for coarse grid
21   USE in_out_manager 
22   USE par_kind
23   USE crslbclnk
24   USE lib_mpp
25   USE ieee_arithmetic   
26
27   IMPLICIT NONE
28
29   PRIVATE
30
31   PUBLIC crs_dom_ope
32   PUBLIC crs_dom_e3, crs_dom_sfc, crs_dom_msk, crs_dom_hgr, crs_dom_coordinates
33   PUBLIC crs_dom_facvol, crs_dom_def, crs_dom_bat, crs_dom_top
34
35   INTERFACE crs_dom_ope
36      MODULE PROCEDURE crs_dom_ope_3d, crs_dom_ope_2d
37   END INTERFACE
38
39   REAL(wp),PUBLIC :: r_inf = 1e+7 !cbr 1e+36
40
41   !! Substitutions
42#  include "domzgr_substitute.h90"
43   
44CONTAINS
45
46
47   SUBROUTINE crs_dom_msk
48   !!===================================================================
49   !
50   !
51   !
52   !!===================================================================
53   INTEGER  ::  ji, jj, jk                   ! dummy loop indices
54   INTEGER  ::  ijis,ijie,ijjs,ijje
55   REAL(wp) ::  zmask
56   !!-------------------------------------------------------------------
57     
58   ! Initialize
59   tmask_crs(:,:,:) = 0.0
60   vmask_crs(:,:,:) = 0.0
61   umask_crs(:,:,:) = 0.0
62   fmask_crs(:,:,:) = 0.0
63   !
64   DO jk = 1, jpkm1
65      DO ji = nldi_crs, nlei_crs
66
67         ijis = mis_crs(ji)
68         ijie = mie_crs(ji)
69
70         DO jj = nldj_crs, nlej_crs
71
72            ijjs = mjs_crs(jj)
73            ijje = mje_crs(jj)
74
75            IF( ijis .le. 0   )WRITE(narea+1000-1,*)"BUG ijis ";CALL FLUSH(narea+1000-1)
76            IF( ijie .gt. jpi )WRITE(narea+1000-1,*)"BUG ijie ";CALL FLUSH(narea+1000-1)
77            IF( ijjs .le. 0   )WRITE(narea+1000-1,*)"BUG ijjs ";CALL FLUSH(narea+1000-1)
78            IF( ijje .gt. jpj )WRITE(narea+1000-1,*)"BUG ijje ";CALL FLUSH(narea+1000-1)
79
80            !IF( jj==nlej_crs .OR. ji==nlei_crs )THEN
81            !WRITE(narea+1000-1,*)"MASK A",ji,jj
82            !WRITE(narea+1000-1,*)"MASK A",ji+nimpp_crs-1,jj+njmpp_crs-1
83            !WRITE(narea+1000-1,*)"MASK B",ijis,ijie,ijjs,ijje
84            !WRITE(narea+1000-1,*)"MASK B",ijis+nimpp-1,ijie+nimpp-1,ijjs+njmpp-1,ijje+njmpp-1
85            !WRITE(narea+1000-1,*)"MASK C",SUM( tmask(ijis:ijie,ijjs:ijje,jk) )
86            !ENDIF
87
88            zmask = 0.0
89            zmask = SUM( tmask(ijis:ijie,ijjs:ijje,jk) )
90            IF ( zmask > 0.0 ) tmask_crs(ji,jj,jk) = 1.0
91
92            zmask = 0.0
93            zmask = SUM( vmask(ijis:ijie,ijje     ,jk) )
94            IF ( zmask > 0.0 ) vmask_crs(ji,jj,jk) = 1.0
95
96            zmask = 0.0
97            zmask = SUM( umask(ijie     ,ijjs:ijje,jk) )
98            IF ( zmask > 0.0 ) umask_crs(ji,jj,jk) = 1.0
99
100            fmask_crs(ji,jj,jk) = fmask(ijie,ijje,jk)
101
102         ENDDO
103      ENDDO
104   ENDDO
105
106   CALL crs_lbc_lnk( tmask_crs, 'T', 1.0 )
107   CALL crs_lbc_lnk( vmask_crs, 'V', 1.0 )
108   CALL crs_lbc_lnk( umask_crs, 'U', 1.0 )
109   CALL crs_lbc_lnk( fmask_crs, 'F', 1.0 )
110
111   !
112   END SUBROUTINE crs_dom_msk
113
114   SUBROUTINE crs_dom_coordinates( p_gphi, p_glam, cd_type, p_gphi_crs, p_glam_crs )
115      !!----------------------------------------------------------------
116      !!               *** SUBROUTINE crs_coordinates ***
117      !! ** Purpose :  Determine the coordinates for the coarse grid
118      !!
119      !! ** Method  :  From the parent grid subset, search for the central
120      !!               point.  For an odd-numbered reduction factor,
121      !!               the coordinate will be that of the central T-cell.
122      !!               For an even-numbered reduction factor, of a non-square
123      !!               coarse grid box, the coordinate will be that of
124      !!               the east or north face or more likely.  For a square
125      !!               coarse grid box, the coordinate will be that of
126      !!               the central f-corner.
127      !!
128      !! ** Input   :  p_gphi = parent grid gphi[t|u|v|f]
129      !!               p_glam = parent grid glam[t|u|v|f]
130      !!               cd_type  = grid type (T,U,V,F)
131      !! ** Output  :  p_gphi_crs = coarse grid gphi[t|u|v|f]
132      !!               p_glam_crs = coarse grid glam[t|u|v|f]
133      !!             
134      !! History. 1 Jun.
135      !!----------------------------------------------------------------
136      !! Arguments
137      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(in)  :: p_gphi  ! Parent grid latitude
138      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(in)  :: p_glam  ! Parent grid longitude
139      CHARACTER(len=1),                     INTENT(in)  :: cd_type   ! grid type (T,U,V,F)
140      REAL(wp), DIMENSION(jpi_crs,jpj_crs), INTENT(out) :: p_gphi_crs  ! Coarse grid latitude
141      REAL(wp), DIMENSION(jpi_crs,jpj_crs), INTENT(out) :: p_glam_crs  ! Coarse grid longitude
142
143      !! Local variables
144      INTEGER :: ji, jj, jk                   ! dummy loop indices
145      INTEGER :: iji, ijj
146      INTEGER  :: ir,jr
147      !!----------------------------------------------------------------
148      p_gphi_crs(:,:)=0._wp
149      p_glam_crs(:,:)=0._wp
150
151 
152      SELECT CASE ( cd_type )
153         CASE ( 'T' )
154            DO jj =  nldj_crs, nlej_crs
155               ijj = mjs_crs(jj) + + INT(0.5*nfacty(jj))
156               DO ji = nldi_crs, nlei_crs
157                  iji = mis_crs(ji) + INT(0.5*nfactx(ji))
158                  p_gphi_crs(ji,jj) = p_gphi(iji,ijj)
159                  p_glam_crs(ji,jj) = p_glam(iji,ijj)
160               ENDDO
161            ENDDO
162         CASE ( 'U' )
163            DO jj =  nldj_crs, nlej_crs
164               ijj = mjs_crs(jj) + INT(0.5*nfacty(jj))
165               DO ji = nldi_crs, nlei_crs
166                  iji = mie_crs(ji)
167                  p_gphi_crs(ji,jj) = p_gphi(iji,ijj)
168                  p_glam_crs(ji,jj) = p_glam(iji,ijj)
169 
170               ENDDO
171            ENDDO
172         CASE ( 'V' )
173            DO jj =  nldj_crs, nlej_crs
174               ijj = mje_crs(jj)
175               DO ji = nldi_crs, nlei_crs
176                  iji = mis_crs(ji) + INT(0.5*nfactx(ji))
177                  p_gphi_crs(ji,jj) = p_gphi(iji,ijj)
178                  p_glam_crs(ji,jj) = p_glam(iji,ijj)
179               ENDDO
180            ENDDO
181         CASE ( 'F' )
182            DO jj =  nldj_crs, nlej_crs
183               ijj = mje_crs(jj)
184               DO ji = nldi_crs, nlei_crs
185                  iji = mie_crs(ji)
186                  p_gphi_crs(ji,jj) = p_gphi(iji,ijj)
187                  p_glam_crs(ji,jj) = p_glam(iji,ijj)
188               ENDDO
189            ENDDO
190      END SELECT
191
192      ! Retroactively add back the boundary halo cells.
193      !CALL crs_lbc_lnk( p_gphi_crs, cd_type, 1.0 )
194      !CALL crs_lbc_lnk( p_glam_crs, cd_type, 1.0 )
195      WRITE(narea+1000-1,*)"end1 glam_crs gphi_crs ",p_glam_crs(1,2),p_gphi_crs(1,2)
196      !
197   END SUBROUTINE crs_dom_coordinates
198
199  SUBROUTINE crs_dom_hgr( p_e1, p_e2, cd_type, p_e1_crs, p_e2_crs )
200      !!----------------------------------------------------------------
201      !!               *** SUBROUTINE crs_dom_hgr ***
202      !!
203      !! ** Purpose :  Get coarse grid horizontal scale factors and unmasked fraction
204      !!
205      !! ** Method  :  For grid types T,U,V,Fthe 2D scale factors of
206      !!               the coarse grid are the sum of the east or north faces of the
207      !!               parent grid subset comprising the coarse grid box.     
208      !!               - e1,e2 Scale factors
209      !!                 Valid arguments:
210      !! ** Inputs  : p_e1, p_e2  = parent grid e1 or e2 (t,u,v,f)
211      !!              cd_type     = grid type (T,U,V,F) for scale factors; for velocities (U or V)
212      !! ** Outputs : p_e1_crs, p_e2_crs  = parent grid e1 or e2 (t,u,v,f)
213      !!
214      !! History.     4 Jun.  Write for WGT and scale factors only
215      !!----------------------------------------------------------------
216      !!
217      !!  Arguments
218      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(in)  :: p_e1     ! Parent grid U,V scale factors (e1)
219      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(in)  :: p_e2     ! Parent grid U,V scale factors (e2)
220      CHARACTER(len=1)                    , INTENT(in)  :: cd_type  ! grid type U,V
221
222      REAL(wp), DIMENSION(jpi_crs,jpj_crs), INTENT(out) :: p_e1_crs ! Coarse grid box 2D quantity
223      REAL(wp), DIMENSION(jpi_crs,jpj_crs), INTENT(out) :: p_e2_crs ! Coarse grid box 2D quantity
224
225      !! Local variables
226      INTEGER :: ji, jj, jk     ! dummy loop indices
227      INTEGER :: ijis,ijie,ijjs,ijje
228      INTEGER :: i1, j1
229 
230      !!---------------------------------------------------------------- 
231      p_e1_crs(:,:)=0._wp
232      p_e2_crs(:,:)=0._wp
233
234      ! Initialize     
235
236         DO ji = nldi_crs, nlei_crs
237
238            ijis = mis_crs(ji)
239            ijie = mie_crs(ji)
240
241            DO jj = nldj_crs, nlej_crs
242
243               ijjs = mjs_crs(jj)
244               ijje = mje_crs(jj)
245
246               i1=INT(0.5*nfactx(ji))
247               j1=INT(0.5*nfacty(jj))
248
249               ! Only for a factro 3 coarsening
250               SELECT CASE ( cd_type )
251                   CASE ( 'T' )
252                      p_e1_crs(ji,jj) = REAL(nn_factx,wp)*p_e1(ijis+i1,ijjs+j1)
253                      p_e2_crs(ji,jj) = REAL(nn_facty,wp)*p_e2(ijis+i1,ijjs+j1)
254                   CASE ( 'U' )
255                      p_e1_crs(ji,jj) = REAL(nn_factx,wp)*p_e1(ijis+i1,ijjs+j1) 
256                      p_e2_crs(ji,jj) = REAL(nn_facty,wp)*p_e2(ijie   ,ijjs+j1) 
257
258                   CASE ( 'V' )
259                      p_e1_crs(ji,jj) = REAL(nn_factx,wp)*p_e1(ijis+i1,ijje   ) 
260                      p_e2_crs(ji,jj) = REAL(nn_facty,wp)*p_e2(ijis+i1,ijjs+j1) 
261                   CASE ( 'F' )
262                      p_e1_crs(ji,jj) = REAL(nn_factx,wp)*p_e1(ijis+i1,ijje   ) 
263                      p_e2_crs(ji,jj) = REAL(nn_facty,wp)*p_e2(ijie   ,ijjs+j1) 
264               END SELECT
265            ENDDO
266         ENDDO
267
268
269      CALL crs_lbc_lnk( p_e1_crs, cd_type, 1.0 )
270      CALL crs_lbc_lnk( p_e2_crs, cd_type, 1.0 )
271
272   END SUBROUTINE crs_dom_hgr
273
274
275   SUBROUTINE crs_dom_facvol( p_mask, cd_type, p_e1, p_e2, p_e3, p_fld1_crs, p_fld2_crs )
276      !!----------------------------------------------------------------
277      !!               *** SUBROUTINE crsfun_wgt ***
278      !! ** Purpose :  Three applications.
279      !!               1) SUM. Get coarse grid horizontal scale factors and unmasked fraction
280      !!               2) VOL. Get coarse grid box volumes
281      !!               3) WGT. Weighting multiplier for volume-weighted and/or
282      !!                       area-weighted averages.
283      !!                       Weights (i.e. the denominator) calculated here
284      !!                       to avoid IF-tests and division.
285      !! ** Method  :  1) SUM.  For grid types T,U,V,F (and W) the 2D scale factors of
286      !!               the coarse grid are the sum of the east or north faces of the
287      !!               parent grid subset comprising the coarse grid box. 
288      !!               The fractions of masked:total surface (3D) on the east,
289      !!               north and top faces is, optionally, also output.
290      !!               - Top face area sum
291      !!                 Valid arguments: cd_type, cd_op='W', p_pmask, p_e1, p_e2
292      !!               - Top face ocean surface fraction
293      !!                 Valid arguments: cd_type, cd_op='W', p_pmask, p_e1, p_e2       
294      !!               - e1,e2 Scale factors
295      !!                 Valid arguments:
296      !!               2) VOL.  For grid types W and T, the coarse grid box
297      !!               volumes are output. Also optionally, the fraction of 
298      !!               masked:total volume of the parent grid subset is output (i.e. facvol).
299      !!               3) WGT. Based on the grid type, the denominator is pre-determined here to 
300      !!               perform area- or volume- weighted averages,
301      !!               to avoid IF-tests and divisions.
302      !! ** Inputs  : p_e1, p_e2  = parent grid e1 or e2 (t,u,v,f)
303      !!              p_pmask     = parent grid mask (T,U,V,F)
304      !!              cd_type     = grid type (T,U,V,F) for scale factors; for velocities (U or V)
305      !!              cd_op       = applied operation (SUM, VOL, WGT)
306      !!              p_fse3      = (Optional) parent grid vertical level thickness (fse3u or fse3v)
307      !! ** Outputs : p_cfield2d_1 = (Optional) 2D field on coarse grid
308      !!              p_cfield2d_2 = (Optional) 2D field on coarse grid
309      !!              p_cfield3d_1 = (Optional) 3D field on coarse grid
310      !!              p_cfield3d_2 = (Optional) 3D field on coarse grid
311      !!
312      !! History.     4 Jun.  Write for WGT and scale factors only
313      !!----------------------------------------------------------------
314      !!
315      !!  Arguments
316      CHARACTER(len=1),                         INTENT(in)  :: cd_type  ! grid type U,V
317      REAL(wp), DIMENSION(jpi,jpj,jpk)        , INTENT(in)  :: p_mask  ! Parent grid U,V mask
318      REAL(wp), DIMENSION(jpi,jpj)            , INTENT(in)  :: p_e1     ! Parent grid U,V scale factors (e1)
319      REAL(wp), DIMENSION(jpi,jpj)            , INTENT(in)  :: p_e2     ! Parent grid U,V scale factors (e2)
320      REAL(wp), DIMENSION(jpi,jpj,jpk)        , INTENT(in)  :: p_e3     ! Parent grid vertical level thickness (fse3u, fse3v)
321
322      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(out) :: p_fld1_crs ! Coarse grid box 3D quantity
323      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(out) :: p_fld2_crs ! Coarse grid box 3D quantity
324
325      !! Local variables
326      REAL(wp)                                :: zdAm
327      INTEGER                                 :: ji, jj, jk
328      INTEGER :: ijis,ijie,ijjs,ijje
329
330      REAL(wp), DIMENSION(:,:,:), POINTER     :: zvol, zmask     
331      !!---------------------------------------------------------------- 
332   
333      CALL wrk_alloc( jpi, jpj, jpk, zvol, zmask )
334
335      p_fld1_crs(:,:,:) = 0.0
336      p_fld2_crs(:,:,:) = 0.0
337
338      DO jk = 1, jpk
339         zvol (:,:,jk) = p_e1(:,:) * p_e2(:,:) * p_e3(:,:,jk) 
340         zmask(:,:,jk) = p_mask(:,:,jk) 
341      ENDDO
342
343      DO jk = 1, jpk
344         DO ji = nldi_crs, nlei_crs
345
346            ijis = mis_crs(ji)
347            ijie = mie_crs(ji)
348
349            DO jj = nldj_crs, nlej_crs
350
351               ijjs = mjs_crs(jj)
352               ijje = mje_crs(jj)
353
354               p_fld1_crs(ji,jj,jk) =  SUM( zvol(ijis:ijie,ijjs:ijje,jk) )
355               zdAm                 =  SUM( zvol(ijis:ijie,ijjs:ijje,jk) * zmask(ijis:ijie,ijjs:ijje,jk) )
356               p_fld2_crs(ji,jj,jk) = zdAm / p_fld1_crs(ji,jj,jk) 
357            ENDDO
358         ENDDO
359      ENDDO
360      !                                             !  Retroactively add back the boundary halo cells.
361      CALL crs_lbc_lnk( p_fld1_crs, cd_type, 1.0 ) 
362      CALL crs_lbc_lnk( p_fld2_crs, cd_type, 1.0 ) 
363      !
364      CALL wrk_dealloc( jpi, jpj, jpk, zvol, zmask )
365      !
366   END SUBROUTINE crs_dom_facvol
367
368
369   SUBROUTINE crs_dom_ope_3d( p_fld, cd_op, cd_type, p_mask, p_fld_crs, p_e12, p_e3, p_surf_crs, p_mask_crs, psgn )
370      !!----------------------------------------------------------------
371      !!               *** SUBROUTINE crsfun_UV ***
372      !! ** Purpose :  Average, area-weighted, of U or V on the east and north faces
373      !!
374      !! ** Method  :  The U and V velocities (3D) are determined as the area-weighted averages
375      !!               on the east and north faces, respectively,
376      !!               of the parent grid subset comprising the coarse grid box.
377      !!               In the case of the V and F grid, the last jrow minus 1 is spurious.
378      !! ** Inputs  : p_e1_e2     = parent grid e1 or e2 (t,u,v,f)
379      !!              cd_type     = grid type (T,U,V,F) for scale factors; for velocities (U or V)
380      !!              psgn        = sign change over north fold (See lbclnk.F90)
381      !!              p_pmask     = parent grid mask (T,U,V,F) for scale factors;
382      !!                                       for velocities (U or V)
383      !!              p_fse3      = parent grid vertical level thickness (fse3u or fse3v)
384      !!              p_pfield    = U or V on the parent grid
385      !!              p_surf_crs  = (Optional) Coarse grid weight for averaging
386      !! ** Outputs : p_cfield3d = 3D field on coarse grid
387      !!
388      !! History.  29 May.  completed draft.
389      !!            4 Jun.  Revision for WGT
390      !!            5 Jun.  Streamline for area-weighted average only ; separate scale factor and weights.
391      !!----------------------------------------------------------------
392      !!
393      !!  Arguments
394      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)        :: p_fld   ! T, U, V or W on parent grid
395      CHARACTER(len=*),                         INTENT(in)           :: cd_op    ! Operation SUM, MAX or MIN
396      CHARACTER(len=1),                         INTENT(in)           :: cd_type    ! grid type U,V
397      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)           :: p_mask    ! Parent grid T,U,V mask
398      REAL(wp), DIMENSION(jpi,jpj),             INTENT(in), OPTIONAL :: p_e12    ! Parent grid T,U,V scale factors (e1 or e2)
399      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in), OPTIONAL :: p_e3     ! Parent grid vertical level thickness (fse3u, fse3v)
400      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(in), OPTIONAL :: p_surf_crs ! Coarse grid area-weighting denominator   
401      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(in), OPTIONAL :: p_mask_crs    ! Coarse grid T,U,V maska
402      REAL(wp),                                 INTENT(in)           :: psgn    ! sign
403
404      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(out)          :: p_fld_crs ! Coarse grid box 3D quantity
405
406      !! Local variables
407      INTEGER  :: ji, jj, jk 
408      INTEGER  :: ijis, ijie, ijjs, ijje
409      INTEGER  :: ini, inj
410      REAL(wp) :: zflcrs, zsfcrs   
411      REAL(wp), DIMENSION(:,:,:), POINTER :: zsurf, zsurfmsk, zmask,ztabtmp
412      INTEGER  :: ir,jr
413      REAL(wp), DIMENSION(nn_factx,nn_facty):: ztmp
414      REAL(wp), DIMENSION(nn_factx*nn_facty):: ztmp1
415      REAL(wp), DIMENSION(:), ALLOCATABLE   :: ztmp2
416      INTEGER , DIMENSION(1)  :: zdim1
417      REAL(wp) :: zmin,zmax
418      !!---------------------------------------------------------------- 
419   
420      p_fld_crs(:,:,:) = 0.0
421
422      SELECT CASE ( cd_op )
423 
424         CASE ( 'VOL' )
425     
426            CALL wrk_alloc( jpi, jpj, jpk, zsurf, zsurfmsk )
427         
428            SELECT CASE ( cd_type )
429           
430               CASE( 'T', 'W','U','V' )
431                  DO jk = 1, jpk
432                     zsurf   (:,:,jk) =  p_e12(:,:) * p_e3(:,:,jk) *  p_mask(:,:,jk) 
433                     zsurfmsk(:,:,jk) =  zsurf(:,:,jk) 
434                  ENDDO
435                  !
436                  DO jk = 1, jpk         
437                     DO jj  = nldj_crs,nlej_crs
438                        ijjs = mjs_crs(jj)
439                        ijje = mje_crs(jj)
440                        DO ji = nldi_crs, nlei_crs
441
442                           ijis = mis_crs(ji)
443                           ijie = mie_crs(ji)
444
445                           zflcrs = SUM( p_fld(ijis:ijie,ijjs:ijje,jk) * zsurfmsk(ijis:ijie,ijjs:ijje,jk) )
446                           zsfcrs = SUM(                                 zsurfmsk(ijis:ijie,ijjs:ijje,jk) )
447
448                           p_fld_crs(ji,jj,jk) = zflcrs
449                           IF( zsfcrs /= 0.0 )  p_fld_crs(ji,jj,jk) = zflcrs / zsfcrs
450                        ENDDO     
451                     ENDDO
452                  ENDDO 
453                  !
454               CASE DEFAULT
455                    STOP
456            END SELECT
457
458            CALL wrk_dealloc( jpi, jpj, jpk, zsurf, zsurfmsk )
459
460         CASE ( 'LOGVOL' )
461
462            CALL wrk_alloc( jpi, jpj, jpk, zsurf, zsurfmsk, ztabtmp )
463
464            ztabtmp(:,:,:)=0._wp
465            WHERE(p_fld* p_mask .NE. 0._wp ) ztabtmp =  LOG10(p_fld * p_mask)*p_mask
466            ztabtmp = ztabtmp * p_mask
467
468            SELECT CASE ( cd_type )
469
470               CASE( 'T', 'W' )
471
472                  DO jk = 1, jpk
473                     zsurf   (:,:,jk) =  p_e12(:,:) * p_e3(:,:,jk) *  p_mask(:,:,jk)
474                     zsurfmsk(:,:,jk) =  zsurf(:,:,jk)
475                  ENDDO
476                  !
477                  DO jk = 1, jpk
478                     DO jj  = nldj_crs,nlej_crs
479                        ijjs = mjs_crs(jj)
480                        ijje = mje_crs(jj)
481                        DO ji = nldi_crs, nlei_crs
482                           ijis = mis_crs(ji)
483                           ijie = mie_crs(ji)
484                           zflcrs = SUM( ztabtmp(ijis:ijie,ijjs:ijje,jk) * zsurfmsk(ijis:ijie,ijjs:ijje,jk) )
485                           zsfcrs = SUM(                                   zsurfmsk(ijis:ijie,ijjs:ijje,jk) )
486                           p_fld_crs(ji,jj,jk) = zflcrs
487                           IF( zsfcrs /= 0.0 )  p_fld_crs(ji,jj,jk) = zflcrs / zsfcrs
488                           p_fld_crs(ji,jj,jk) = 10 ** ( p_fld_crs(ji,jj,jk) *  p_mask_crs(ji,jj,jk) ) * p_mask_crs(ji,jj,jk)
489                        ENDDO
490                     ENDDO
491                  ENDDO
492               CASE DEFAULT
493                    STOP
494            END SELECT
495
496            CALL wrk_dealloc( jpi, jpj, jpk, zsurf, zsurfmsk ,ztabtmp )
497
498         CASE ( 'MED' )
499
500            CALL wrk_alloc( jpi, jpj, jpk, zsurf, zsurfmsk )
501
502            SELECT CASE ( cd_type )
503
504               CASE( 'T', 'W' )
505                  DO jk = 1, jpk
506                     zsurf   (:,:,jk) =  p_e12(:,:) * p_e3(:,:,jk) *  p_mask(:,:,jk)
507                     zsurfmsk(:,:,jk) =  zsurf(:,:,jk)
508                  ENDDO
509                  !
510                  DO jk = 1, jpk
511
512                     DO jj  = nldj_crs,nlej_crs
513
514                        ijjs = mjs_crs(jj)
515                        ijje = mje_crs(jj)
516                        inj  = ijje-ijjs+1
517
518                        DO ji = nldi_crs, nlei_crs
519                           ijis = mis_crs(ji)
520                           ijie = mie_crs(ji)
521                           ini  = ijie-ijis+1
522
523                           ztmp(1:ini,1:inj)= p_fld(ijis:ijie,ijjs:ijje,jk)
524                           zdim1(1) = nn_factx*nn_facty
525                           ztmp1(:) = RESHAPE( ztmp(:,:) , zdim1 )
526                           CALL PIKSRT(nn_factx*nn_facty,ztmp1)
527
528                           ir=0
529                           jr=1
530                           DO WHILE( jr .LE. nn_factx*nn_facty )
531                              IF( ztmp1(jr) == 0. ) THEN
532                                 ir=jr
533                                 jr=jr+1
534                              ELSE
535                                 EXIT
536                              ENDIF
537                           ENDDO
538                           IF( ir .LE. nn_factx*nn_facty-1 )THEN
539                              ALLOCATE( ztmp2(nn_factx*nn_facty-ir) )
540                              ztmp2(1:nn_factx*nn_facty-ir) = ztmp1(1+ir:nn_factx*nn_facty)
541                              jr=INT( 0.5 * REAL(nn_factx*nn_facty-ir,wp) )+1
542                              p_fld_crs(ji,jj,jk) = ztmp2(jr)
543                              DEALLOCATE( ztmp2 )
544                           ELSE
545                              p_fld_crs(ji,jj,jk) = 0._wp
546                           ENDIF
547
548                        ENDDO
549                     ENDDO
550                  ENDDO
551               CASE DEFAULT
552                    STOP
553            END SELECT
554
555           CALL wrk_dealloc( jpi, jpj, jpk, zsurf, zsurfmsk )
556 
557         CASE ( 'SUM' )
558         
559            CALL wrk_alloc( jpi, jpj, jpk, zsurfmsk )
560
561            IF( PRESENT( p_e3 ) ) THEN
562               DO jk = 1, jpk
563                  zsurfmsk(:,:,jk) =  p_e12(:,:) * p_e3(:,:,jk) * p_mask(:,:,jk) 
564               ENDDO
565            ELSE
566               DO jk = 1, jpk
567                  zsurfmsk(:,:,jk) =  p_e12(:,:) * p_mask(:,:,jk) 
568               ENDDO
569            ENDIF
570
571            SELECT CASE ( cd_type )
572           
573               CASE( 'T', 'W' )
574       
575                  DO jk = 1, jpk
576                     DO jj  = nldj_crs,nlej_crs
577                        ijjs = mjs_crs(jj)
578                        ijje = mje_crs(jj)
579                        DO ji = nldi_crs, nlei_crs
580                           ijis = mis_crs(ji)
581                           ijie = mie_crs(ji)
582
583                           p_fld_crs(ji,jj,jk) = SUM( p_fld(ijis:ijie,ijjs:ijje,jk) * zsurfmsk(ijis:ijie,ijjs:ijje,jk) )
584                        ENDDO
585                     ENDDO
586                  ENDDO
587
588               CASE( 'V' )
589
590
591                  DO jk = 1, jpk
592                     DO jj  = nldj_crs,nlej_crs
593                        ijjs = mjs_crs(jj)
594                        ijje = mje_crs(jj)
595                        DO ji = nldi_crs, nlei_crs
596                           ijis = mis_crs(ji)
597                           ijie = mie_crs(ji)
598
599                           p_fld_crs(ji,jj,jk) = SUM( p_fld(ijis:ijie,ijje,jk) * zsurfmsk(ijis:ijie,ijje,jk) )
600                        ENDDO
601                     ENDDO
602                  ENDDO
603
604               CASE( 'U' )
605
606                  DO jk = 1, jpk
607                     DO jj  = nldj_crs,nlej_crs
608                        ijjs = mjs_crs(jj)
609                        ijje = mje_crs(jj)
610                        DO ji = nldi_crs, nlei_crs
611                           ijis = mis_crs(ji)
612                           ijie = mie_crs(ji)
613
614                           p_fld_crs(ji,jj,jk) = SUM( p_fld(ijie,ijjs:ijje,jk) * zsurfmsk(ijie,ijjs:ijje,jk) )
615                        ENDDO
616                     ENDDO
617                  ENDDO
618
619              END SELECT
620
621              IF( PRESENT( p_surf_crs ) ) THEN
622                 WHERE ( p_surf_crs /= 0.0 ) p_fld_crs(:,:,:) = p_fld_crs(:,:,:) / p_surf_crs(:,:,:)
623              ENDIF
624
625              CALL wrk_dealloc( jpi, jpj, jpk, zsurfmsk )
626
627         CASE ( 'MAX' )    !  search the max of unmasked grid cells
628
629            CALL wrk_alloc( jpi, jpj, jpk, zmask )
630
631            DO jk = 1, jpk
632               zmask(:,:,jk) = p_mask(:,:,jk) 
633            ENDDO
634
635            SELECT CASE ( cd_type )
636           
637               CASE( 'T', 'W' )
638       
639                  DO jk = 1, jpk
640                     DO jj  = nldj_crs,nlej_crs
641                        ijjs = mjs_crs(jj)
642                        ijje = mje_crs(jj)
643                        DO ji = nldi_crs, nlei_crs
644                           ijis = mis_crs(ji)
645                           ijie = mie_crs(ji)
646                           p_fld_crs(ji,jj,jk) = MAXVAL( p_fld(ijis:ijie,ijjs:ijje,jk) * zmask(ijis:ijie,ijjs:ijje,jk) - &
647                                                       & ( ( 1._wp - zmask(ijis:ijie,ijjs:ijje,jk))* r_inf )                )
648                        ENDDO
649                     ENDDO
650                  ENDDO
651 
652               CASE( 'V' )
653                  CALL ctl_stop('MAX operator and V case not available')
654           
655               CASE( 'U' )
656                  CALL ctl_stop('MAX operator and U case not available')
657
658            END SELECT
659
660            CALL wrk_dealloc( jpi, jpj, jpk, zmask )
661
662         CASE ( 'MIN' )      !   Search the min of unmasked grid cells
663
664            CALL wrk_alloc( jpi, jpj, jpk, zmask )
665            DO jk = 1, jpk
666               zmask(:,:,jk) = p_mask(:,:,jk)
667            ENDDO
668
669            SELECT CASE ( cd_type )
670
671               CASE( 'T', 'W' )
672
673                  DO jk = 1, jpk
674                     DO jj  = nldj_crs,nlej_crs
675                        ijjs = mjs_crs(jj)
676                        ijje = mje_crs(jj)
677                        DO ji = nldi_crs, nlei_crs
678                           ijis = mis_crs(ji)
679                           ijie = mie_crs(ji)
680
681                           p_fld_crs(ji,jj,jk) = MINVAL( p_fld(ijis:ijie,ijjs:ijje,jk) * zmask(ijis:ijie,ijjs:ijje,jk) + &
682                                                       & ( 1._wp - zmask(ijis:ijie,ijjs:ijje,jk)* r_inf )                )
683                        ENDDO
684                     ENDDO
685                  ENDDO
686
687           
688               CASE( 'V' )
689                  CALL ctl_stop('MIN operator and V case not available')
690           
691               CASE( 'U' )
692                  CALL ctl_stop('MIN operator and U case not available')
693         
694            END SELECT
695            !
696            CALL wrk_dealloc( jpi, jpj, jpk, zmask )
697            !
698         END SELECT
699         !
700         CALL crs_lbc_lnk( p_fld_crs, cd_type, psgn )
701         !
702    END SUBROUTINE crs_dom_ope_3d
703
704    SUBROUTINE crs_dom_ope_2d( p_fld, cd_op, cd_type, p_mask, p_fld_crs, p_e12, p_e3, p_surf_crs, p_mask_crs, psgn )
705      !!----------------------------------------------------------------
706      !!               *** SUBROUTINE crsfun_UV ***
707      !! ** Purpose :  Average, area-weighted, of U or V on the east and north faces
708      !!
709      !! ** Method  :  The U and V velocities (3D) are determined as the area-weighted averages
710      !!               on the east and north faces, respectively,
711      !!               of the parent grid subset comprising the coarse grid box.
712      !!               In the case of the V and F grid, the last jrow minus 1 is spurious.
713      !! ** Inputs  : p_e1_e2     = parent grid e1 or e2 (t,u,v,f)
714      !!              cd_type     = grid type (T,U,V,F) for scale factors; for velocities (U or V)
715      !!              psgn        = sign change over north fold (See lbclnk.F90)
716      !!              p_pmask     = parent grid mask (T,U,V,F) for scale factors;
717      !!                                       for velocities (U or V)
718      !!              p_fse3      = parent grid vertical level thickness (fse3u or fse3v)
719      !!              p_pfield    = U or V on the parent grid
720      !!              p_surf_crs  = (Optional) Coarse grid weight for averaging
721      !! ** Outputs : p_cfield3d = 3D field on coarse grid
722      !!
723      !! History.  29 May.  completed draft.
724      !!            4 Jun.  Revision for WGT
725      !!            5 Jun.  Streamline for area-weighted average only ; separate scale factor and weights.
726      !!----------------------------------------------------------------
727      !!
728      !!  Arguments
729      REAL(wp), DIMENSION(jpi,jpj),             INTENT(in)           :: p_fld   ! T, U, V or W on parent grid
730      CHARACTER(len=*),                         INTENT(in)           :: cd_op    ! Operation SUM, MAX or MIN
731      CHARACTER(len=1),                         INTENT(in)           :: cd_type    ! grid type U,V
732      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)           :: p_mask    ! Parent grid T,U,V mask
733      REAL(wp), DIMENSION(jpi,jpj),             INTENT(in), OPTIONAL :: p_e12    ! Parent grid T,U,V scale factors (e1 or e2)
734      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in), OPTIONAL :: p_e3     ! Parent grid vertical level thickness (fse3u, fse3v)
735      REAL(wp), DIMENSION(jpi_crs,jpj_crs)    , INTENT(in), OPTIONAL :: p_surf_crs ! Coarse grid area-weighting denominator   
736      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(in), OPTIONAL :: p_mask_crs    ! Coarse grid T,U,V mask
737      REAL(wp),                                 INTENT(in)           :: psgn   
738
739      REAL(wp), DIMENSION(jpi_crs,jpj_crs)    , INTENT(out)          :: p_fld_crs ! Coarse grid box 3D quantity
740
741      !! Local variables
742      INTEGER  :: ji, jj, jk                 ! dummy loop indices
743      INTEGER ::  ijis, ijie, ijjs, ijje
744      REAL(wp) :: zflcrs, zsfcrs   
745      REAL(wp), DIMENSION(:,:), POINTER :: zsurfmsk   
746
747      !!---------------------------------------------------------------- 
748 
749      p_fld_crs(:,:) = 0.0
750
751      SELECT CASE ( TRIM(cd_op) )
752     
753        CASE ( 'VOL' )
754
755            CALL wrk_alloc( jpi, jpj, zsurfmsk )
756            zsurfmsk(:,:) =  p_e12(:,:) * p_e3(:,:,1) * p_mask(:,:,1)
757
758            DO jj  = nldj_crs,nlej_crs
759               ijjs = mjs_crs(jj)
760               ijje = mje_crs(jj)
761               DO ji = nldi_crs, nlei_crs
762                  ijis = mis_crs(ji)
763                  ijie = mie_crs(ji)
764
765                  zflcrs = SUM( p_fld(ijis:ijie,ijjs:ijje) * zsurfmsk(ijis:ijie,ijjs:ijje) )
766                  zsfcrs = SUM(                              zsurfmsk(ijis:ijie,ijjs:ijje) )
767
768                  p_fld_crs(ji,jj) = zflcrs
769                  IF( zsfcrs /= 0.0 )  p_fld_crs(ji,jj) = zflcrs / zsfcrs
770               ENDDO
771            ENDDO
772            CALL wrk_dealloc( jpi, jpj, zsurfmsk )
773            !
774
775         CASE ( 'SUM' )
776         
777            CALL wrk_alloc( jpi, jpj, zsurfmsk )
778            IF( PRESENT( p_e3 ) ) THEN
779               zsurfmsk(:,:) =  p_e12(:,:) * p_e3(:,:,1) * p_mask(:,:,1)
780            ELSE
781               zsurfmsk(:,:) =  p_e12(:,:) * p_mask(:,:,1)
782            ENDIF
783
784            SELECT CASE ( cd_type )
785
786               CASE( 'T', 'W' )
787
788                  DO jj  = nldj_crs,nlej_crs
789                     ijjs = mjs_crs(jj)
790                     ijje = mje_crs(jj)
791                     DO ji = nldi_crs, nlei_crs
792                        ijis = mis_crs(ji)
793                        ijie = mie_crs(ji)
794                        p_fld_crs(ji,jj) = SUM( p_fld(ijis:ijie,ijjs:ijje) * zsurfmsk(ijis:ijie,ijjs:ijje) )
795                     ENDDO
796                  ENDDO
797           
798               CASE( 'V' )
799
800                  DO jj  = nldj_crs,nlej_crs
801                     ijjs = mjs_crs(jj)
802                     ijje = mje_crs(jj)
803                     DO ji = nldi_crs, nlei_crs
804                        ijis = mis_crs(ji)
805                        ijie = mie_crs(ji)
806                        p_fld_crs(ji,jj) = SUM( p_fld(ijis:ijie,ijje) * zsurfmsk(ijis:ijie,ijje) )
807                     ENDDO
808                  ENDDO
809
810               CASE( 'U' )
811
812                  DO jj  = nldj_crs,nlej_crs
813                     ijjs = mjs_crs(jj)
814                     ijje = mje_crs(jj)
815                     DO ji = nldi_crs, nlei_crs
816                        ijis = mis_crs(ji)
817                        ijie = mie_crs(ji)
818                        p_fld_crs(ji,jj) = SUM( p_fld(ijie,ijjs:ijje) * zsurfmsk(ijie,ijjs:ijje) )
819                     ENDDO
820                  ENDDO
821
822              END SELECT
823
824              IF( PRESENT( p_surf_crs ) ) THEN
825                 WHERE ( p_surf_crs /= 0.0 ) p_fld_crs(:,:) = p_fld_crs(:,:) / p_surf_crs(:,:)
826              ENDIF
827
828              CALL wrk_dealloc( jpi, jpj, zsurfmsk )
829
830         CASE ( 'MAX' )
831
832            SELECT CASE ( cd_type )
833           
834               CASE( 'T', 'W' )
835 
836                  DO jj  = nldj_crs,nlej_crs
837                     ijjs = mjs_crs(jj)
838                     ijje = mje_crs(jj)
839                     DO ji = nldi_crs, nlei_crs
840                        ijis = mis_crs(ji)
841                        ijie = mie_crs(ji)
842                        p_fld_crs(ji,jj) = MAXVAL( p_fld(ijis:ijie,ijjs:ijje) * p_mask(ijis:ijie,ijjs:ijje,1) - &
843                                                 & ( 1._wp - p_mask(ijis:ijie,ijjs:ijje,1) )                    )
844                     ENDDO
845                  ENDDO
846           
847               CASE( 'V' )
848                  CALL ctl_stop('MAX operator and V case not available')
849           
850               CASE( 'U' )
851                  CALL ctl_stop('MAX operator and U case not available')
852
853              END SELECT
854
855         CASE ( 'MIN' )      !   Search the min of unmasked grid cells
856
857           SELECT CASE ( cd_type )
858
859              CASE( 'T', 'W' )
860
861                  DO jj  = nldj_crs,nlej_crs
862                     ijjs = mjs_crs(jj)
863                     ijje = mje_crs(jj)
864                     DO ji = nldi_crs, nlei_crs
865                        ijis = mis_crs(ji)
866                        ijie = mie_crs(ji)
867                        p_fld_crs(ji,jj) = MINVAL( p_fld(ijis:ijie,ijjs:ijje) * p_mask(ijis:ijie,ijjs:ijje,1) + &
868                                                 & ( 1._wp - p_mask(ijis:ijie,ijjs:ijje,1) )                    )
869                     ENDDO
870                  ENDDO
871           
872               CASE( 'V' )
873                  CALL ctl_stop('MIN operator and V case not available')
874           
875               CASE( 'U' )
876                  CALL ctl_stop('MIN operator and U case not available')
877
878              END SELECT
879             !
880         END SELECT
881         !
882         CALL crs_lbc_lnk( p_fld_crs, cd_type, psgn )
883         !
884   END SUBROUTINE crs_dom_ope_2d
885
886   SUBROUTINE crs_dom_e3( p_e1, p_e2, p_e3, p_sfc_2d_crs,  p_sfc_3d_crs, cd_type, p_mask, p_e3_crs, p_e3_max_crs)
887      !!---------------------------------------------------------------- 
888      !!
889      !!
890      !!
891      !!
892      !!----------------------------------------------------------------
893      !!  Arguments
894      CHARACTER(len=1),                         INTENT(in)          :: cd_type           ! grid type T, W ( U, V, F)
895      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)          :: p_mask            ! Parent grid T mask
896      REAL(wp), DIMENSION(jpi,jpj)    ,         INTENT(in)          :: p_e1, p_e2        ! 2D tracer T or W on parent grid
897      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)          :: p_e3              ! 3D tracer T or W on parent grid
898      REAL(wp), DIMENSION(jpi_crs,jpj_crs)    , INTENT(in),OPTIONAL :: p_sfc_2d_crs      ! Coarse grid box east or north face quantity
899      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(in),OPTIONAL :: p_sfc_3d_crs      ! Coarse grid box east or north face quantity
900      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(inout)       :: p_e3_crs          ! Coarse grid box east or north face quantity
901      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(inout)       :: p_e3_max_crs      ! Coarse grid box east or north face quantity
902
903      !! Local variables
904      INTEGER ::  ji, jj, jk                   ! dummy loop indices
905      INTEGER ::  ijis, ijie, ijjs, ijje 
906      REAL(wp) :: ze3crs 
907
908      !!---------------------------------------------------------------- 
909      p_e3_crs    (:,:,:) = 0._wp
910      p_e3_max_crs(:,:,:) = 0._wp
911   
912
913      SELECT CASE ( cd_type )
914
915         CASE ('T')
916
917            DO jk = 1, jpk
918               DO ji = nldi_crs, nlei_crs
919
920                  ijis = mis_crs(ji)
921                  ijie = mie_crs(ji)
922
923                  DO jj = nldj_crs, nlej_crs
924
925                     ijjs = mjs_crs(jj)
926                     ijje = mje_crs(jj)
927
928                     p_e3_max_crs(ji,jj,jk) = MAXVAL( p_e3(ijis:ijie,ijjs:ijje,jk) * p_mask(ijis:ijie,ijjs:ijje,jk) )
929
930                     ze3crs = SUM( p_e1(ijis:ijie,ijjs:ijje) * p_e2(ijis:ijie,ijjs:ijje) * p_e3(ijis:ijie,ijjs:ijje,jk) * p_mask(ijis:ijie,ijjs:ijje,jk) )
931                     IF( p_sfc_3d_crs(ji,jj,jk) .NE. 0._wp )p_e3_crs(ji,jj,jk) = ze3crs / p_sfc_3d_crs(ji,jj,jk)
932
933                  ENDDO
934               ENDDO
935            ENDDO
936
937         CASE ('U')
938
939            DO jk = 1, jpk
940               DO ji = nldi_crs, nlei_crs
941
942                  ijis = mis_crs(ji)
943                  ijie = mie_crs(ji)
944
945                  DO jj = nldj_crs, nlej_crs
946
947                     ijjs = mjs_crs(jj)
948                     ijje = mje_crs(jj)
949
950                     p_e3_max_crs(ji,jj,jk) = MAXVAL( p_e3(ijie,ijjs:ijje,jk) * p_mask(ijie,ijjs:ijje,jk) )
951
952                     ze3crs = SUM( p_e2(ijie,ijjs:ijje) * p_e3(ijie,ijjs:ijje,jk) * p_mask(ijie,ijjs:ijje,jk) )
953                     IF( p_sfc_2d_crs(ji,jj) .NE. 0._wp )p_e3_crs(ji,jj,jk) = ze3crs / p_sfc_2d_crs(ji,jj)
954                  ENDDO
955               ENDDO
956            ENDDO
957
958         CASE ('V')
959
960            DO jk = 1, jpk
961               DO ji = nldi_crs, nlei_crs
962
963                  ijis = mis_crs(ji)
964                  ijie = mie_crs(ji)
965
966                  DO jj = nldj_crs, nlej_crs
967
968                     ijjs = mjs_crs(jj)
969                     ijje = mje_crs(jj)
970
971                     p_e3_max_crs(ji,jj,jk) = MAXVAL( p_e3(ijis:ijie,ijje,jk) * p_mask(ijis:ijie,ijje,jk) )
972
973                     ze3crs = SUM( p_e1(ijis:ijie,ijje) * p_e3(ijis:ijie,ijje,jk) * p_mask(ijis:ijie,ijje,jk) )
974                     IF( p_sfc_2d_crs(ji,jj) .NE. 0._wp )p_e3_crs(ji,jj,jk) = ze3crs / p_sfc_2d_crs(ji,jj)
975
976                  ENDDO
977               ENDDO
978            ENDDO
979
980         CASE ('W')
981
982            DO jk = 1, jpk
983               DO ji = nldi_crs, nlei_crs
984
985                  ijis = mis_crs(ji)
986                  ijie = mie_crs(ji)
987
988                  DO jj = nldj_crs, nlej_crs
989
990                     ijjs = mjs_crs(jj)
991                     ijje = mje_crs(jj)
992
993                     p_e3_max_crs(ji,jj,jk) = MAXVAL( p_e3(ijis:ijie,ijjs:ijje,jk) * p_mask(ijis:ijie,ijjs:ijje,jk) )
994
995                     ze3crs = SUM( p_e1(ijis:ijie,ijjs:ijje) * p_e2(ijis:ijie,ijjs:ijje) * p_e3(ijis:ijie,ijjs:ijje,jk) * p_mask(ijis:ijie,ijjs:ijje,jk) )
996                     IF( p_sfc_3d_crs(ji,jj,jk) .NE. 0._wp )p_e3_crs(ji,jj,jk) = ze3crs / p_sfc_3d_crs(ji,jj,jk)
997
998                  ENDDO
999               ENDDO
1000            ENDDO
1001
1002      END SELECT
1003
1004      CALL crs_lbc_lnk( p_e3_max_crs, cd_type, 1.0 )
1005      CALL crs_lbc_lnk( p_e3_crs    , cd_type, 1.0 )
1006
1007   END SUBROUTINE crs_dom_e3
1008
1009   SUBROUTINE crs_dom_sfc(p_mask, cd_type, p_surf_crs, p_surf_crs_msk,  p_e1, p_e2, p_e3 )
1010      !!=========================================================================================
1011      !!
1012      !!
1013      !!=========================================================================================
1014      !!  Arguments
1015      CHARACTER(len=1),                         INTENT(in)           :: cd_type      ! grid type T, W ( U, V, F)
1016      REAL(wp), DIMENSION(jpi,jpj,jpk)        , INTENT(in)           :: p_mask       ! Parent grid T mask
1017      REAL(wp), DIMENSION(jpi,jpj)            , INTENT(in), OPTIONAL :: p_e1, p_e2         ! 3D tracer T or W on parent grid
1018      REAL(wp), DIMENSION(jpi,jpj,jpk)        , INTENT(in), OPTIONAL :: p_e3         ! 3D tracer T or W on parent grid
1019      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(out)          :: p_surf_crs ! Coarse grid box east or north face quantity
1020      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT(out)          :: p_surf_crs_msk ! Coarse grid box east or north face quantity
1021
1022      !! Local variables
1023      INTEGER  :: ji, jj, jk                   ! dummy loop indices
1024      INTEGER  :: ijis,ijie,ijjs,ijje
1025      REAL(wp), DIMENSION(:,:,:), POINTER :: zsurf, zsurfmsk   
1026      !!---------------------------------------------------------------- 
1027      ! Initialize
1028      p_surf_crs(:,:,:)=0._wp
1029      p_surf_crs_msk(:,:,:)=0._wp
1030
1031      CALL wrk_alloc( jpi, jpj, jpk, zsurf, zsurfmsk )
1032      !
1033      SELECT CASE ( cd_type )
1034     
1035         CASE ('W')   
1036            DO jk = 1, jpk
1037               zsurf(:,:,jk) = p_e1(:,:) * p_e2(:,:) 
1038            ENDDO
1039
1040         CASE ('V')     
1041            DO jk = 1, jpk
1042               zsurf(:,:,jk) = p_e1(:,:) * p_e3(:,:,jk) 
1043            ENDDO
1044 
1045         CASE ('U')     
1046            DO jk = 1, jpk
1047               zsurf(:,:,jk) = p_e2(:,:) * p_e3(:,:,jk) 
1048            ENDDO
1049
1050         CASE DEFAULT
1051            DO jk = 1, jpk
1052               zsurf(:,:,jk) = p_e1(:,:) * p_e2(:,:) 
1053            ENDDO
1054      END SELECT
1055
1056      DO jk = 1, jpk
1057         zsurfmsk(:,:,jk) = zsurf(:,:,jk) * p_mask(:,:,jk)
1058      ENDDO
1059
1060      SELECT CASE ( cd_type )
1061
1062         CASE ('W')
1063
1064            DO jk = 1, jpk
1065               DO jj = nldj_crs,nlej_crs
1066                  ijjs=mjs_crs(jj)
1067                  ijje=mje_crs(jj)
1068                  DO ji = nldi_crs,nlei_crs
1069                     ijis=mis_crs(ji)
1070                     ijie=mie_crs(ji)
1071                     p_surf_crs    (ji,jj,jk) =  SUM(zsurf   (ijis:ijie,ijjs:ijje,jk) )
1072                     p_surf_crs_msk(ji,jj,jk) =  SUM(zsurfmsk(ijis:ijie,ijjs:ijje,jk) )
1073                  ENDDO     
1074               ENDDO
1075            ENDDO   
1076
1077         CASE ('U')
1078
1079            DO jk = 1, jpk
1080               DO jj = nldj_crs,nlej_crs
1081                  ijjs=mjs_crs(jj)
1082                  ijje=mje_crs(jj)
1083                  DO ji = nldi_crs,nlei_crs
1084                     ijis=mis_crs(ji)
1085                     ijie=mie_crs(ji)
1086                     p_surf_crs    (ji,jj,jk) =  SUM(zsurf   (ijie,ijjs:ijje,jk) )
1087                     p_surf_crs_msk(ji,jj,jk) =  SUM(zsurfmsk(ijie,ijjs:ijje,jk) )
1088                  ENDDO
1089               ENDDO
1090            ENDDO
1091
1092         CASE ('V')
1093
1094            DO jk = 1, jpk
1095               DO jj = nldj_crs,nlej_crs
1096                  ijjs=mjs_crs(jj)
1097                  ijje=mje_crs(jj)
1098                  DO ji = nldi_crs,nlei_crs
1099                     ijis=mis_crs(ji)
1100                     ijie=mie_crs(ji)
1101                     p_surf_crs    (ji,jj,jk) =  SUM(zsurf   (ijis:ijie,ijje,jk) )
1102                     p_surf_crs_msk(ji,jj,jk) =  SUM(zsurfmsk(ijis:ijie,ijje,jk) )
1103                  ENDDO
1104               ENDDO
1105            ENDDO
1106
1107      END SELECT
1108
1109      CALL crs_lbc_lnk( p_surf_crs    , cd_type , 1. )
1110      CALL crs_lbc_lnk( p_surf_crs_msk, cd_type , 1. )
1111
1112      CALL wrk_dealloc( jpi, jpj, jpk, zsurf, zsurfmsk )
1113
1114   END SUBROUTINE crs_dom_sfc
1115   
1116   SUBROUTINE crs_dom_def
1117      !!----------------------------------------------------------------
1118      !!               *** SUBROUTINE crs_dom_def ***
1119      !! ** Purpose :  Three applications.
1120      !!               1) Define global domain indice of the croasening grid
1121      !!               2) Define local domain indice of the croasening grid
1122      !!               3) Define the processor domain indice for a croasening grid
1123      !!----------------------------------------------------------------
1124      INTEGER :: ji,jj,jk,ijjgloT,ijis,ijie,ijjs,ijje,jn      ! dummy indices
1125      INTEGER :: ierr                                ! allocation error status
1126      INTEGER :: iproci,iprocj,iproc,iprocno,iprocso,iimppt_crs
1127      INTEGER :: ii_start,ii_end,ij_start,ij_end
1128      !!----------------------------------------------------------------
1129 
1130 
1131      !==============================================================================================
1132      ! Define global and local domain sizes
1133      !==============================================================================================
1134      SELECT CASE ( jperio )
1135
1136      CASE ( 0, 1 )
1137
1138      !exp for a regional configuration
1139      ! i=1 => ghost cell 
1140
1141              !--!--!--!--!--!--!--!--!--!
1142              !        !        !        !
1143              !        !        !        !
1144              !        !        !        !
1145              !--!--!--!--!--!--!--!--!--!
1146      !IBI36               1  2  3  4  5 
1147      !IBI12  !    1   !   2    !   3    !   
1148 
1149      jpiglo_crs   = INT( (jpiglo - 4) / nn_factx ) + 4
1150      jpjglo_crs   = INT( (jpjglo - 4) / nn_facty ) + 4
1151
1152      jpj_crs = ( jpjglo_crs -1  - 2 * jprecj + (jpnj-1) ) / jpnj + 2 * jprecj
1153      IF( MOD( jpjglo_crs   - 2 * jprecj + (jpnj-1) , jpnj ) .NE.  0._wp )THEN
1154         jpj_crs = jpj_crs + 1
1155      ENDIF
1156
1157      CASE ( 3, 4 )   
1158      jpiglo_crs   = INT( (jpiglo - 2) / nn_factx ) + 2
1159      jpjglo_crs   = INT( (jpjglo - MOD(jpjglo, nn_facty)) / nn_facty ) + 2
1160
1161      jpj_crs = ( jpjglo_crs -1  - 2 * jprecj + (jpnj-1) ) / jpnj + 2 * jprecj
1162
1163      END SELECT
1164     
1165      jpiglo_crsm1 = jpiglo_crs - 1
1166      jpjglo_crsm1 = jpjglo_crs - 1 
1167
1168      jpi_crs = ( jpiglo_crs   - 2 * jpreci + (jpni-1) ) / jpni + 2 * jpreci
1169       
1170      jpi_crsm1   = jpi_crs - 1
1171      jpj_crsm1   = jpj_crs - 1
1172
1173      npolj_crs   = npolj
1174     
1175      ierr = crs_dom_alloc()          ! allocate most coarse grid arrays
1176
1177      !==============================================================================================
1178      ! Define processor domain indices
1179      !==============================================================================================
1180      IF( .NOT. lk_mpp ) THEN
1181
1182         nimpp_crs  = 1
1183         njmpp_crs  = 1
1184         nlci_crs   = jpi_crs
1185         nlcj_crs   = jpj_crs
1186         nldi_crs   = 1
1187         nldj_crs   = 1
1188         nlei_crs   = jpi_crs
1189         nlej_crs   = jpj_crs
1190
1191      ELSE
1192
1193         nimpp_crs  = 1
1194         njmpp_crs  = 1
1195         nlci_crs   = jpi_crs
1196         nlcj_crs   = jpj_crs
1197         nldi_crs   = 1
1198         nldj_crs   = 1
1199         nlei_crs   = jpi_crs
1200         nlej_crs   = jpj_crs
1201
1202         !==============================================================================================
1203         ! mpp_ini2
1204         !==============================================================================================
1205
1206         !order of local domain in i and j directions
1207         DO ji = 1 , jpni
1208            DO jj = 1 ,jpnj
1209               IF( nfipproc(ji,jj)  == narea-1 )THEN
1210                  iproci=ji
1211                  iprocj=jj
1212               ENDIF
1213            ENDDO
1214         ENDDO
1215
1216         !==========================================================================
1217         ! check
1218         !==========================================================================
1219         !WRITE(narea+1000-1,*)"nfipproc(ji,jj),narea :",nfipproc(iproci,iprocj),narea
1220         !WRITE(narea+1000-1,*)"proc i,j ",iproci,iprocj
1221         !WRITE(narea+1000-1,*)"jpni  jpnj jpnij ",jpni,jpnj,jpnij
1222         !WRITE(narea+1000-1,*)"nperio jperio ",nperio,jperio
1223         !WRITE(narea+1000-1,*)"nowe noea",nowe,noea
1224         !WRITE(narea+1000-1,*)"noso nono",noso,nono
1225         !WRITE(narea+1000-1,*)"nbondi nbondj ",nbondi,nbondj
1226         !WRITE(narea+1000-1,*)"jpiglo jpjglo ",jpiglo,jpjglo
1227         !WRITE(narea+1000-1,*)"jpi jpj ",jpi,jpj
1228         !WRITE(narea+1000-1,*)"nbondi nbondj",nbondi,nbondj
1229         !WRITE(narea+1000-1,*)"nimpp njmpp ",nimpp,njmpp
1230         !WRITE(narea+1000-1,*)"loc jpi nldi,nlei,nlci ",jpi, nldi        ,nlei         ,nlci
1231         !WRITE(narea+1000-1,*)"glo jpi nldi,nlei      ",jpi, nldi+nimpp-1,nlei+nimpp-1
1232         !WRITE(narea+1000-1,*)"loc jpj nldj,nlej,nlcj ",jpj, nldj        ,nlej         ,nlcj
1233         !WRITE(narea+1000-1,*)"glo jpj nldj,nlej      ",jpj, nldj+njmpp-1,nlej+njmpp-1
1234         !WRITE(narea+1000-1,*)"jpiglo_crs jpjglo_crs ",jpiglo_crs,jpjglo_crs
1235         !WRITE(narea+1000-1,*)"jpi_crs jpj_crs ",jpi_crs,jpj_crs
1236         !WRITE(narea+1000-1,*)"glamt gphit ",glamt(1,1),gphit(jpi,jpj),glamt(jpi,jpj),gphit(jpi,jpj)
1237         !WRITE(narea+1000-1,*)"min max tmask ",MINVAL(tmask),MAXVAL(tmask)
1238         !CALL FLUSH(narea+1000-1)
1239         !==========================================================================
1240         ! coarsened domain: dimensions along I
1241         !==========================================================================
1242
1243         !------------------------------------------------------------------------------------
1244         !I-1 fill mis2_crs and mie2_crs: arrays to switch from physic grid to coarsened grid
1245         !------------------------------------------------------------------------------------
1246
1247         ! !--------!--------!--------!
1248         ! !        !        !        !
1249         ! !        !        !        !
1250         ! !        !        !        !
1251         ! !--------!--------!--------!
1252         ! !        !        !        !
1253         ! !        ! ji,jj  !        !
1254         ! !        !        !        !
1255         ! !--------!--------!--------!
1256         ! !        !        !        !
1257         ! !        !        !        !
1258         ! !        !        !        !
1259         ! !--------!--------!--------!
1260         !  mis2_crs(ji)      mie2_crs(ji)
1261       
1262
1263         SELECT CASE ( jperio )
1264
1265         CASE ( 0, 1 )
1266
1267            DO ji=1,jpiglo_crs
1268               ijis=nn_factx*ji-6
1269               ijie=nn_factx*ji-4
1270               mis2_crs(ji)=ijis
1271               mie2_crs(ji)=ijie
1272               !WRITE(narea+1000-1,*)"glo crs",ji,ijis,ijie,ijis-nimpp+1,ijie-nimpp+1
1273            ENDDO
1274
1275         CASE ( 3, 4 )    !   3, 4 : T-Pivot at North Fold: make correspondinf the pivot points of the 2 grids
1276
1277            DO ji=1,jpiglo_crs
1278               ijis=nn_factx*(ji-1)-2
1279               ijie=nn_factx*(ji-1)
1280               mis2_crs(ji)=ijis
1281               mie2_crs(ji)=ijie
1282               !WRITE(narea+1000-1,*)"glo crs",ji,ijis,ijie,ijis-nimpp+1,ijie-nimpp+1
1283            ENDDO
1284
1285         CASE DEFAULT
1286            WRITE(numout,*) 'crs_init. Only jperio = 0, 1, 3, 4 supported; narea: ',narea
1287         END SELECT
1288
1289         !-------------------------------------------------------------------------------
1290         ! I-2 find the first CRS cell which is inside the physic grid inner domain
1291         !-------------------------------------------------------------------------------
1292         ! ijis           : global indice of the first CRS cell which inside the physic grid inner domain
1293         ! mis2_crs(ijis) : global indice of the bottom-left physic cell corresponding to ijis cell
1294         ! ii_start       : local  ndice of the bottom-left physic cell corresponding to ijis cell
1295
1296         ji=1
1297         DO WHILE( mis2_crs(ji) - nimpp + 1 .LT. 1 ) 
1298            ji=ji+1
1299            IF( ji==jpiglo_crs )EXIT
1300         END DO
1301
1302         ijis=ji
1303         ii_start = mis2_crs(ijis)-nimpp+1
1304         !WRITE(narea+1000-1,*)"start ",ijis,mis2_crs(ijis),ii_start ; CALL FLUSH(narea+1000-1)
1305
1306         !----------------------------------------------------------------------------------------------
1307         ! I-3 compute nldi_crs
1308         !     special case  when the domain on the left is land:
1309         !         compute mis2_crs and mie2_crs for the first cell of the local domain
1310         !---------------------------------------------------------------------------------------------
1311         nldi_crs = 2
1312         IF( nowe == -1 .AND. ( jperio==3 .OR. jperio==4 .OR. jperio==0 .OR. jperio==1 ) )THEN
1313
1314            mie2_crs(ijis-1) = mis2_crs(ijis)-1
1315             
1316            SELECT CASE(ii_start)
1317               CASE(1)
1318                  nldi_crs=2
1319                  mie2_crs(ijis-1) = 1
1320                  mis2_crs(ijis-1) = 1
1321               CASE(2)
1322                  nldi_crs=2
1323                  mis2_crs(ijis-1) = mie2_crs(ijis-1)
1324               CASE(3)
1325                  nldi_crs=2
1326                  mis2_crs(ijis-1) = mie2_crs(ijis-1) -1
1327               !CASE DEFAULT
1328               !  WRITE(narea+8000-1,*)"WRONG VALUE FOR iistart ",ii_start
1329            END SELECT
1330
1331         ENDIF
1332         
1333         !----------------------------------------------------------------------------------------------
1334         ! I-4 compute nimpp_crs
1335         !---------------------------------------------------------------------------------------------
1336         nimpp_crs = ijis-1
1337         IF( nimpp==1 )nimpp_crs=1
1338
1339         !WRITE(narea+1000-1,*)" nimpp_crs ",nimpp_crs
1340
1341         !-------------------------------------------------------------------------------
1342         ! I-5 find the last CRS cell which is inside the physic grid inner domain
1343         !-------------------------------------------------------------------------------
1344         ! ijie           : global indice of the last CRS cell which inside the physic grid inner domain
1345
1346         ji=jpiglo_crs
1347         DO WHILE( mie2_crs(ji) - nimpp + 1 .GT. nlci )
1348            ji=ji-1
1349            IF( ji==1 )EXIT
1350         END DO
1351         ijie=ji
1352         !WRITE(narea+1000-1,*)"end ",ijie ; CALL FLUSH(narea+1000-1)
1353
1354         !-------------------------------------------------------------------------------
1355         ! I-6 compute nlei_crs and nlci_crs
1356         !-------------------------------------------------------------------------------
1357         nlei_crs = ijie-nimpp_crs+1
1358         nlci_crs = nlei_crs+1
1359
1360         !-------------------------------------------------------------------------------
1361         ! I-7 local to global and global to local indices for CRS grid
1362         !-------------------------------------------------------------------------------
1363         DO ji = 1, jpi_crs
1364            mig_crs(ji) = ji + nimpp_crs - 1
1365            !WRITE(narea+1000-1,*)"crs loctoglo",ji,mig_crs(ji) ; CALL FLUSH(narea+1000-1)
1366         ENDDO
1367         DO ji = 1, jpiglo_crs
1368            mi0_crs(ji) = MAX( 1, MIN( ji - nimpp_crs + 1 , jpi_crs + 1 ) )
1369            mi1_crs(ji) = MAX( 0, MIN( ji - nimpp_crs + 1 , jpi_crs     ) )
1370         ENDDO
1371
1372         !---------------------------------------------------------------------------------------
1373         ! I-8 CRS to physic grid: local indices mis_crs and mie_crs and local coarsening factor
1374         !---------------------------------------------------------------------------------------
1375         DO ji = 1, nlei_crs
1376            mis_crs(ji) = mis2_crs(mig_crs(ji)) - nimpp + 1
1377            mie_crs(ji) = mie2_crs(mig_crs(ji)) - nimpp + 1
1378            nfactx(ji)  = mie_crs(ji)-mis_crs(ji)+1
1379         ENDDO
1380
1381         !-----------------------------------------------------------------------------
1382         !adapt western indices for E-W BC ( jperio 3/4 ) or ghost cells ( jperio 0/1)
1383         !-----------------------------------------------------------------------------
1384         IF( iproci == 1 ) THEN
1385
1386            SELECT CASE ( jperio )
1387
1388            CASE ( 3, 4 )
1389
1390               nldi_crs=1
1391               mis_crs(1)  = 1
1392               mie_crs(1)  = 1
1393               mis2_crs(1) = 1
1394               mie2_crs(1) = 1
1395
1396            CASE ( 0, 1 )
1397               nldi_crs=1
1398               !the first cell is a ghsot cell; so set indices to 1 just to have an positive indice
1399               mis_crs(1)  = 1
1400               mie_crs(1)  = 1
1401               mis2_crs(1) = 1
1402               mie2_crs(1) = 1
1403               nfactx(1)   = 1
1404
1405               !the second coarsened cell correspond to the first (ghost) and second cell of HR grid
1406               mis_crs(2)  = 1
1407               mie_crs(2)  = 2
1408               mis2_crs(2) = 1
1409               mie2_crs(2) = 2
1410               nfactx(2)   = 2
1411
1412            CASE DEFAULT
1413               WRITE(numout,*) 'crs_init. Only jperio = 0, 1, 3, 4 supported; narea: ',narea
1414            END SELECT
1415       
1416         ENDIF
1417
1418         !----------------------------------------------------------------------------
1419         !adapt eastern indices for E-W BC ( jperio 3/4 ) or ghost cells ( jperio 0/1)
1420         !----------------------------------------------------------------------------
1421         SELECT CASE ( jperio )
1422
1423         CASE ( 3, 4 )
1424
1425            IF( iproci == jpni ) THEN
1426               nlei_crs           = jpiglo_crs-nimpp_crs+1
1427               nlci_crs           = nlei_crs
1428               mis_crs(nlei_crs)  = 1
1429               mie_crs(nlei_crs)  = 1
1430               mis2_crs(nlei_crs) = 1
1431               mie2_crs(nlei_crs) = 1
1432               nfactx(nlei_crs)   = 0
1433            ENDIF
1434
1435         CASE ( 0, 1 )
1436            IF( iproci == jpni ) THEN
1437
1438             !              |                 |                 |                 |
1439             !            __|_____!_____!_____|_____!_____!_____|_____!_____!_____|
1440             !HR    grid loc                         nlei
1441             !HR    grid glo                        jpiglo
1442             !HRCRS grid loc                      nlei_crs-1           nlei
1443             !HRCRS grid glo                      jpiglo-1            jpiglo
1444
1445
1446               !last cell ( = ghost cell)
1447               nlei_crs           = jpiglo_crs-nimpp_crs+1
1448               nlci_crs           = nlei_crs
1449     
1450               mis2_crs(nlei_crs) = jpiglo
1451               mie2_crs(nlei_crs) = jpiglo
1452               mis_crs (nlei_crs) = jpiglo     - nimpp + 1
1453               mie_crs (nlei_crs) = jpiglo     - nimpp + 1
1454               nfactx  (nlei_crs) = 1
1455
1456               mis2_crs(nlei_crs-1) = jpiglo - 1
1457               mie2_crs(nlei_crs-1) = jpiglo
1458               mis_crs (nlei_crs-1) = jpiglo - 1 - nimpp + 1 
1459               mie_crs (nlei_crs-1) = jpiglo     - nimpp + 1
1460               nfactx  (nlei_crs-1) = 2
1461
1462            ENDIF
1463         CASE DEFAULT
1464            WRITE(numout,*) 'crs_init. Only jperio = 0, 1, 3, 4 supported; narea: ',narea
1465         END SELECT
1466
1467         
1468         !WRITE(narea+1000-1,*)"loc crs jpi nldi,nlei,nlci ",jpi_crs, nldi_crs            ,nlei_crs             ,nlci_crs
1469         !CALL FLUSH(narea+1000-1)
1470         !WRITE(narea+1000-1,*)"glo crs jpi nldi,nlei      ",jpi_crs, nldi_crs+nimpp_crs-1,nlei_crs+nimpp_crs-1
1471         !CALL FLUSH(narea+1000-1)
1472         !DO ji = 1, jpi_crs
1473         !   WRITE(narea+1000-1,'(A15,7(1X,I4.4))')"crs test i",ji,ji+nimpp_crs-1,mis_crs(ji),mie_crs(ji), &
1474         !        &                                        mis_crs(ji)+nimpp-1,mie_crs(ji)+nimpp-1,nfactx(ji)
1475         !ENDDO
1476
1477         !==========================================================================
1478         ! coarsened domain: dimensions along J
1479         !==========================================================================
1480
1481         !-----------------------------------------------------------------------------------
1482         !J-1 fill mjs2_crs and mje2_crs: arrays to switch from physic grid to coarsened grid
1483         !-----------------------------------------------------------------------------------
1484
1485         ! !--------!--------!--------!
1486         ! !        !        !        !
1487         ! !        !        !        !
1488         ! !        !        !        ! mje2_crs(jj)
1489         ! !--------!--------!--------!
1490         ! !        !        !        !
1491         ! !        ! ji,jj  !        !
1492         ! !        !        !        !
1493         ! !--------!--------!--------!
1494         ! !        !        !        !
1495         ! !        !        !        ! mjs2_crs(jj)
1496         ! !        !        !        !
1497         ! !--------!--------!--------!
1498
1499         SELECT CASE ( jperio )
1500
1501         CASE ( 0, 1 )    !
1502
1503            DO jj=1,jpjglo_crs
1504               ijjs=nn_facty*jj-6
1505               ijje=nn_facty*jj-4
1506               mjs2_crs(jj)=ijjs
1507               mje2_crs(jj)=ijje
1508               !WRITE(narea+1000-1,*)"glo crs",jj,ijjs,ijje,ijjs-njmpp+1,ijje-njmpp+1
1509            ENDDO
1510
1511         CASE ( 3, 4 )    !   3, 4 : T-Pivot at North Fold
1512
1513            DO jj=1,jpjglo_crs
1514               ijjs=nn_facty*(jj) - 5 + MOD(jpjglo, nn_facty)
1515               ijje=nn_facty*(jj) - 3 + MOD(jpjglo, nn_facty)
1516               mjs2_crs(jj)=ijjs
1517               mje2_crs(jj)=ijje
1518               !WRITE(narea+1000-1,*)"glo crs",jj,ijjs,ijje,ijjs-njmpp+1,ijje-njmpp+1
1519            ENDDO
1520
1521         CASE DEFAULT
1522            WRITE(numout,*) 'crs_init. Only jperio = 0, 1, 3, 4 supported; narea: ',narea
1523         END SELECT
1524
1525         !-------------------------------------------------------------------------------
1526         ! J-2 find the first CRS cell which is inside the physic grid inner domain
1527         !-------------------------------------------------------------------------------
1528         ! ijjs           : global indice of the first CRS cell which inside the physic grid inner domain
1529         ! mis2_crs(ijjs) : global indice of the bottom-left physic cell corresponding to ijis cell
1530         ! ij_start       : local  ndice of the bottom-left physic cell corresponding to ijis cell
1531
1532         jj=1
1533         DO WHILE( mjs2_crs(jj) - njmpp + 1 .LT. 1 )
1534            jj=jj+1
1535            IF( jj==jpjglo_crs )EXIT
1536         END DO
1537
1538         ijjs=jj
1539         ij_start = mjs2_crs(ijjs)-njmpp+1
1540
1541         !----------------------------------------------------------------------------------------------
1542         ! J-3 compute nldj_crs
1543         !     global: special case when the domain on the north is land:
1544         !             compute mjs2_crs and mje2_crs for the first cell of the local domain
1545         !     local :  special case for southern domain
1546         !---------------------------------------------------------------------------------------------
1547         nldj_crs = 2
1548
1549         IF( jperio==3 .OR. jperio==4 )THEN
1550
1551            IF( noso == -1 )THEN
1552
1553               mje2_crs(ijjs-1) = mjs2_crs(ijjs)-1
1554
1555               SELECT CASE(ij_start)
1556                  CASE(1)
1557                     nldj_crs=2
1558                     mje2_crs(ijjs-1) = -1
1559                     mjs2_crs(ijjs-1) = -1
1560                  CASE(2)
1561                     nldj_crs=2
1562                     mjs2_crs(ijjs-1) = mje2_crs(ijjs-1)
1563                  CASE(3)
1564                     nldj_crs=2
1565                     mjs2_crs(ijjs-1) = mje2_crs(ijjs-1) -1
1566                  !CASE DEFAULT
1567                  !   WRITE(narea+8000-1,*)"WRONG VALUE FOR iistart ",ii_start
1568               END SELECT
1569
1570            ENDIF
1571         ENDIF
1572         IF( jperio==0 .OR. jperio==1 )THEN
1573            IF( njmpp==1 )nldj_crs=1
1574         END IF
1575
1576         !----------------------------------------------------------------------------------------------
1577         ! J-4 compute nimpp_crs
1578         !---------------------------------------------------------------------------------------------
1579         njmpp_crs = ijjs-1
1580         IF( njmpp==1 )njmpp_crs=1
1581         !WRITE(narea+1000-1,*)" njmpp_crs ",njmpp_crs
1582
1583         !-------------------------------------------------------------------------------
1584         ! J-5 find the last CRS cell which is inside the physic grid inner domain
1585         !-------------------------------------------------------------------------------
1586         ! ijje           : global indice of the last CRS cell which inside the physic grid inner domain
1587
1588         jj=jpjglo_crs
1589         DO WHILE( mje2_crs(jj) - njmpp + 1 .GT. nlcj )
1590            jj=jj-1
1591            IF( jj==1 )EXIT
1592         END DO
1593         ijje=jj
1594
1595         !-------------------------------------------------------------------------------
1596         ! J-6 compute nlej_crs and nlcj_crs
1597         !-------------------------------------------------------------------------------
1598         nlej_crs=ijje-njmpp_crs+1
1599         nlcj_crs=nlej_crs+jprecj
1600
1601         !north-fold BC for global
1602         IF( iprocj == jpnj )THEN
1603            IF( jperio==3 .OR. jperio==4 )THEN
1604               nlej_crs=nlej_crs+1
1605               nlcj_crs=nlej_crs
1606            ENDIF
1607         ENDIF
1608
1609         !-------------------------------------------------------------------------------
1610         ! J-7 local to global and global to local indices for CRS grid
1611         !-------------------------------------------------------------------------------
1612         DO jj = 1, jpj_crs
1613            mjg_crs(jj) = jj + njmpp_crs - 1
1614         ENDDO
1615         DO jj = 1, jpjglo_crs
1616            mj0_crs(jj) = MAX( 1, MIN( jj - njmpp_crs + 1 , jpj_crs + 1 ) )
1617            mj1_crs(jj) = MAX( 0, MIN( jj - njmpp_crs + 1 , jpj_crs     ) )
1618         ENDDO
1619
1620         !---------------------------------------------------------------------------------------
1621         ! J-8 CRS to physic grid: local indices mis_crs and mie_crs and local coarsening factor
1622         !---------------------------------------------------------------------------------------
1623         DO jj = 1, nlej_crs
1624            mjs_crs(jj) = mjs2_crs(mjg_crs(jj)) - njmpp + 1
1625            mje_crs(jj) = mje2_crs(mjg_crs(jj)) - njmpp + 1
1626            nfacty(jj)  = mje_crs(jj)-mjs_crs(jj)+1
1627            !WRITE(narea+1000-1,*)"test J",jj,mjg_crs(jj),mjs_crs(jj),mje_crs(jj),mjs_crs(jj)+njmpp-1,mje_crs(jj)+njmpp-1,nfacty(jj)
1628         ENDDO
1629
1630         !-----------------------------------------------------------------------------
1631         !adpat indices for southern side of the domain (ghost cells for jperio 0/1 )
1632         !-----------------------------------------------------------------------------
1633         IF( iprocj == 1 ) THEN
1634
1635            SELECT CASE ( jperio )
1636
1637            CASE ( 3, 4 )
1638               ! do nothing
1639            CASE ( 0, 1 )
1640                nldj_crs = 1
1641               !the first cell is a ghsot cell; so set indices to 1 just to have an positive indice
1642               mjs_crs(1)  = 1
1643               mje_crs(1)  = 1
1644               mjs2_crs(1) = 1
1645               mje2_crs(1) = 1
1646               nfacty(1)   = 1
1647
1648               !the second coarsened cell correspond to the first (ghost) and second cell of HR grid
1649               mjs_crs(2)  = 1
1650               mje_crs(2)  = 2
1651               mjs2_crs(2) = 1
1652               mje2_crs(2) = 2
1653               nfacty(2)   = 2
1654
1655            CASE DEFAULT
1656               WRITE(numout,*) 'crs_init. Only jperio = 0, 1, 3, 4 supported; narea: ',narea
1657            END SELECT
1658
1659         ENDIF
1660
1661         !-----------------------------------------------
1662         !adapt northen indices for ghost cells ( jperio 0/1)
1663         !-----------------------------------------------
1664         SELECT CASE ( jperio )
1665
1666         CASE ( 3, 4 )
1667            IF( iprocj == jpnj  .AND. jj == nlej_crs )THEN
1668               mjs_crs(jj) = nlej
1669               mjs2_crs(mjg_crs(jj)) = mjg(nlej)
1670               mje_crs(jj) = nlej
1671               mje2_crs(mjg_crs(jj)) = mjg(nlej)
1672            ENDIF
1673
1674         CASE ( 0, 1 )
1675            IF( iprocj == jpnj ) THEN
1676               !north side of the domain
1677
1678               !last cell ( = ghost cell)
1679               nlej_crs           = jpjglo_crs-njmpp_crs+1
1680               nlcj_crs           = nlej_crs
1681
1682               mjs2_crs(nlej_crs) = jpjglo
1683               mje2_crs(nlej_crs) = jpjglo
1684               mjs_crs (nlej_crs) = jpjglo     - njmpp + 1
1685               mje_crs (nlej_crs) = jpjglo     - njmpp + 1
1686               nfacty  (nlej_crs) = 1
1687
1688               mjs2_crs(nlej_crs-1) = jpjglo - 1
1689               mje2_crs(nlej_crs-1) = jpjglo
1690               mjs_crs (nlej_crs-1) = jpjglo - 1 - njmpp + 1
1691               mje_crs (nlej_crs-1) = jpjglo     - njmpp + 1
1692               nfacty  (nlej_crs-1) = 2
1693
1694            ENDIF
1695         CASE DEFAULT
1696           WRITE(numout,*) 'crs_init. Only jperio = 0, 1, 3, 4 supported; narea: ',narea
1697         END SELECT
1698
1699         !WRITE(narea+1000-1,*)"loc crs jpj nldj,nlej,nlcj ",jpj_crs, nldj_crs ,nlej_crs,nlcj_crs ; CALL FLUSH(narea+1000-1)
1700         !WRITE(narea+1000-1,*)"glo crs jpj nldj,nlej      ",jpj_crs, nldj_crs+njmpp_crs-1,nlej_crs+njmpp_crs-1 ;CALL FLUSH(narea+1000-1)
1701         !DO jj = 1, jpj_crs
1702         !   WRITE(narea+1000-1,'(A15,7(1X,I4.4))')"crs test j",jj,jj+njmpp_crs-1,mjs_crs(jj),mje_crs(jj), &
1703         !        & mjs_crs(jj)+njmpp-1,mje_crs(jj)+njmpp-1,nfacty(jj)
1704         !ENDDO
1705
1706         !==========================================================================
1707         ! send local start and end indices to all procs
1708         !==========================================================================
1709
1710         nldit_crs(:)=0 ; nleit_crs(:)=0 ; nlcit_crs(:)=0 ; nimppt_crs(:)=0
1711         nldjt_crs(:)=0 ; nlejt_crs(:)=0 ; nlcjt_crs(:)=0 ; njmppt_crs(:)=0
1712
1713         CALL mppgatheri((/nlci_crs/),0,nlcit_crs) ; CALL mppgatheri((/nlcj_crs/),0,nlcjt_crs) 
1714         CALL mppgatheri((/nldi_crs/),0,nldit_crs) ; CALL mppgatheri((/nldj_crs/),0,nldjt_crs) 
1715         CALL mppgatheri((/nlei_crs/),0,nleit_crs) ; CALL mppgatheri((/nlej_crs/),0,nlejt_crs) 
1716         CALL mppgatheri((/nimpp_crs/),0,nimppt_crs) ; CALL mppgatheri((/njmpp_crs/),0,njmppt_crs) 
1717
1718         DO jj = 1 ,jpnj
1719            DO ji = 1 , jpni
1720               jn=nfipproc(ji,jj)+1
1721               IF( jn .GE. 1 )THEN
1722                  nfiimpp_crs(ji,jj)=nimppt_crs(jn)
1723               ELSE
1724                  nfiimpp_crs(ji,jj) = ANINT( REAL( (nfiimpp(ji,jj) + 1 ) / nn_factx, wp ) ) + 1
1725               ENDIF
1726            ENDDO
1727         ENDDO
1728 
1729         !nogather=T
1730         nfsloop_crs = 1
1731         nfeloop_crs = nlci_crs
1732         DO jn = 2,jpni-1
1733            IF( nfipproc(jn,jpnj) .eq. (narea - 1) )THEN
1734               IF (nfipproc(jn - 1 ,jpnj) .eq. -1 )THEN
1735                  nfsloop_crs = nldi_crs
1736               ENDIF
1737               IF( nfipproc(jn + 1,jpnj) .eq. -1 )THEN
1738                  nfeloop_crs = nlei_crs
1739               ENDIF
1740            ENDIF
1741         END DO
1742
1743         !==========================================================================
1744         ! check
1745         !==========================================================================
1746         !WRITE(narea+1000-1,*)"mpp_crs ",nimpp_crs,njmpp_crs
1747         !WRITE(narea+1000-1,*)"loc crs jpi nldi,nlei,nlci ",jpi_crs, nldi_crs            ,nlei_crs             ,nlci_crs
1748         !WRITE(narea+1000-1,*)"glo crs jpi nldi,nlei      ",jpi_crs, nldi_crs+nimpp_crs-1,nlei_crs+nimpp_crs-1
1749         !WRITE(narea+1000-1,*)"loc crs jpj nldj,nlej,nlcj ",jpj_crs, nldj_crs            ,nlej_crs             ,nlcj_crs
1750         !WRITE(narea+1000-1,*)"glo crs jpj nldj,nlej      ",jpj_crs, nldj_crs+njmpp_crs-1,nlej_crs+njmpp_crs-1
1751
1752         !==========================================================================
1753         ! Save the parent grid information
1754         !==========================================================================
1755         IF( jpizoom /= 1 .OR. jpjzoom /= 1)    STOP  !cbr mettre un ctlstp et ailleurs ( crsini )
1756         jpi_full    = jpi
1757         jpj_full    = jpj
1758         jpim1_full  = jpim1
1759         jpjm1_full  = jpjm1
1760         npolj_full  = npolj
1761         jpiglo_full = jpiglo
1762         jpjglo_full = jpjglo
1763
1764         nlcj_full   = nlcj
1765         nlci_full   = nlci
1766         nldi_full   = nldi
1767         nldj_full   = nldj
1768         nlei_full   = nlei
1769         nlej_full   = nlej
1770         nimpp_full  = nimpp     
1771         njmpp_full  = njmpp
1772     
1773         nlcit_full(:)  = nlcit(:)
1774         nldit_full(:)  = nldit(:)
1775         nleit_full(:)  = nleit(:)
1776         nimppt_full(:) = nimppt(:)
1777         nlcjt_full(:)  = nlcjt(:)
1778         nldjt_full(:)  = nldjt(:)
1779         nlejt_full(:)  = nlejt(:)
1780         njmppt_full(:) = njmppt(:)
1781     
1782         nfsloop_full = nfsloop
1783         nfeloop_full = nfeloop
1784
1785         nfiimpp_full(:,:) = nfiimpp(:,:) 
1786
1787
1788         !==========================================================================
1789         ! control
1790         !==========================================================================
1791         CALL dom_grid_crs  !swich from mother grid to coarsened grid
1792
1793         IF(lwp) THEN
1794            WRITE(numout,*)
1795            WRITE(numout,*) 'crs_init : coarse grid dimensions'
1796            WRITE(numout,*) '~~~~~~~   coarse domain global j-dimension           jpjglo = ', jpjglo
1797            WRITE(numout,*) '~~~~~~~   coarse domain global i-dimension           jpiglo = ', jpiglo
1798            WRITE(numout,*) '~~~~~~~   coarse domain local  i-dimension              jpi = ', jpi
1799            WRITE(numout,*) '~~~~~~~   coarse domain local  j-dimension              jpj = ', jpj
1800            WRITE(numout,*)
1801            WRITE(numout,*) ' nproc  = '     , nproc
1802            WRITE(numout,*) ' nlci   = '     , nlci
1803            WRITE(numout,*) ' nlcj   = '     , nlcj
1804            WRITE(numout,*) ' nldi   = '     , nldi
1805            WRITE(numout,*) ' nldj   = '     , nldj
1806            WRITE(numout,*) ' nlei   = '     , nlei
1807            WRITE(numout,*) ' nlej   = '     , nlej
1808            WRITE(numout,*) ' nlei_full='    , nlei_full
1809            WRITE(numout,*) ' nldi_full='    , nldi_full
1810            WRITE(numout,*) ' nimpp  = '     , nimpp
1811            WRITE(numout,*) ' njmpp  = '     , njmpp
1812            WRITE(numout,*) ' njmpp_full  = ', njmpp_full
1813            WRITE(numout,*)
1814         ENDIF
1815     
1816         CALL dom_grid_glo ! switch from coarsened grid to mother grid
1817     
1818         nrestx = MOD( nn_factx, 2 )   ! check if even- or odd- numbered reduction factor
1819         nresty = MOD( nn_facty, 2 )
1820
1821         IF( nresty == 0 )THEN
1822            IF( npolj == 3 ) npolj_crs = 5
1823            IF( npolj == 5 ) npolj_crs = 3
1824         ENDIF     
1825     
1826         rfactxy = nn_factx * nn_facty
1827     
1828      ENDIF ! lk_mpp
1829      !
1830      nistr = mis_crs(2)  ;   niend = mis_crs(nlci_crs - 1)
1831      njstr = mjs_crs(3)  ;   njend = mjs_crs(nlcj_crs - 1)
1832      !
1833      !
1834   END SUBROUTINE crs_dom_def
1835   
1836   SUBROUTINE crs_dom_bat
1837      !!----------------------------------------------------------------
1838      !!               *** SUBROUTINE crs_dom_bat ***
1839      !! ** Purpose :  coarsenig bathy
1840      !!----------------------------------------------------------------
1841      !!
1842      !!  local variables
1843      INTEGER  :: ji,jj,jk      ! dummy indices
1844      REAL(wp), DIMENSION(:,:)  , POINTER :: zmbk
1845      !!----------------------------------------------------------------
1846   
1847      CALL wrk_alloc( jpi_crs, jpj_crs, zmbk )
1848   
1849      mbathy_crs(:,:) = jpkm1
1850      mbkt_crs(:,:) = 1
1851      mbku_crs(:,:) = 1
1852      mbkv_crs(:,:) = 1
1853
1854
1855      DO jj = 1, jpj_crs
1856         DO ji = 1, jpi_crs
1857            jk = 0
1858            DO WHILE( tmask_crs(ji,jj,jk+1) > 0.)
1859               jk = jk + 1
1860            ENDDO
1861            mbathy_crs(ji,jj) = float( jk )
1862         ENDDO
1863      ENDDO
1864     
1865      CALL wrk_alloc( jpi_crs, jpj_crs, zmbk )
1866
1867      zmbk(:,:) = 0.0
1868      zmbk(:,:) = REAL( mbathy_crs(:,:), wp ) ;   CALL crs_lbc_lnk(zmbk,'T',1.0)   ;   mbathy_crs(:,:) = INT( zmbk(:,:) )
1869
1870
1871      !
1872      IF(lwp) WRITE(numout,*)
1873      IF(lwp) WRITE(numout,*) '    crsini : mbkt is ocean bottom k-index of T-, U-, V- and W-levels '
1874      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
1875      !
1876      mbkt_crs(:,:) = MAX( mbathy_crs(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
1877      !                                     ! bottom k-index of W-level = mbkt+1
1878
1879      DO jj = 1, jpj_crsm1                      ! bottom k-index of u- (v-) level
1880         DO ji = 1, jpi_crsm1
1881            mbku_crs(ji,jj) = MIN(  mbkt_crs(ji+1,jj  ) , mbkt_crs(ji,jj)  )
1882            mbkv_crs(ji,jj) = MIN(  mbkt_crs(ji  ,jj+1) , mbkt_crs(ji,jj)  )
1883         END DO
1884      END DO
1885
1886      ! convert into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
1887      zmbk(:,:) = 1.e0;   
1888      zmbk(:,:) = REAL( mbku_crs(:,:), wp )   ;   CALL crs_lbc_lnk(zmbk,'U',1.0) ; mbku_crs  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
1889      zmbk(:,:) = REAL( mbkv_crs(:,:), wp )   ;   CALL crs_lbc_lnk(zmbk,'V',1.0) ; mbkv_crs  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
1890      !
1891      CALL wrk_dealloc( jpi_crs, jpj_crs, zmbk )
1892      !
1893   END SUBROUTINE crs_dom_bat
1894
1895   SUBROUTINE crs_dom_top
1896      !!----------------------------------------------------------------
1897      !!               *** SUBROUTINE crs_dom_top ***
1898      !!
1899      !! ** Purpose :  coarsening mikt
1900      !!
1901      !!----------------------------------------------------------------
1902      INTEGER ::  ji, jj, jk
1903      INTEGER ::  ijis, ijie, ijjs, ijje
1904
1905      !!----------------------------------------------------------------
1906      mikt_crs(:,:)=1
1907
1908      DO jj  = nldj_crs,nlej_crs
1909
1910         ijjs = mjs_crs(jj)
1911         ijje = mje_crs(jj)
1912
1913         DO ji = nldi_crs, nlei_crs
1914
1915            ijis = mis_crs(ji)
1916            ijie = mie_crs(ji)
1917
1918            mikt_crs(ji,jj) = MINVAL( mikt(ijis:ijie,ijjs:ijje) )
1919         ENDDO
1920      ENDDO
1921
1922      mikt_crs = MAX( mikt_crs , 1 )
1923
1924   END SUBROUTINE crs_dom_top
1925
1926   SUBROUTINE PIKSRT(N,ARR)
1927     INTEGER                  ,INTENT(IN) :: N
1928     REAL(kind=8),DIMENSION(N),INTENT(INOUT) :: ARR
1929
1930     INTEGER      :: i,j
1931     REAL(kind=8) :: a
1932     !!----------------------------------------------------------------
1933
1934     DO j=2, N
1935       a=ARR(j)
1936       DO i=j-1,1,-1
1937          IF(ARR(i)<=a) goto 10
1938          ARR(i+1)=ARR(i)
1939       ENDDO
1940       i=0
194110     ARR(i+1)=a
1942     ENDDO
1943     RETURN
1944
1945   END SUBROUTINE PIKSRT
1946
1947
1948END MODULE crsdom
Note: See TracBrowser for help on using the repository browser.