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

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

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

Last change on this file since 10210 was 10210, checked in by cbricaud, 6 years ago

modifications for limited area configuration

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