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 branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS – NEMO

source: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdom.F90 @ 8096

Last change on this file since 8096 was 8096, checked in by cbricaud, 7 years ago

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