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

Last change on this file since 10171 was 10171, checked in by cbricaud, 22 months ago

3.6/coarsening branch: init mikt_crs to avoid segmentation fault when initialise BCG tracer with netcdf file

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