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.
crsini.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/crsini.F90 @ 10107

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

add call to mikt coarsening in nemo3.6 coarsening branch
  • Property svn:keywords set to Id
File size: 15.6 KB
RevLine 
[4015]1MODULE crsini   
2   !!======================================================================
3   !!                         ***  MODULE crsini  ***
4   !!            Manage the grid coarsening module initialization
5   !!======================================================================
6   !!  History     2012-05   (J. Simeon, G. Madec, C. Ethe, C. Calone) Original code
7   !!----------------------------------------------------------------------
8
9   USE timing                   ! Timing
10   USE par_oce                  ! For parameter jpi,jpj,jphgr_msh
11   USE dom_oce                  ! For parameters in par_oce (jperio, lk_vvl)
12   USE crs                  ! Coarse grid domain
13   USE phycst, ONLY: omega, rad ! physical constants
14   USE wrk_nemo 
15   USE in_out_manager
16   USE par_kind, ONLY: wp
17   USE iom
18   USE crsdom
19   USE crsdomwri
20   USE crslbclnk
21   USE lib_mpp
[5601]22   USE ldftra_crs
[6772]23   USE ieee_arithmetic
[4015]24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC  crs_init
29
30   !! * Substitutions
31#  include "domzgr_substitute.h90"
32
33CONTAINS
34   
35   SUBROUTINE crs_init 
36      !!-------------------------------------------------------------------
37      !!                     *** SUBROUTINE crs_init
38      !!  ** Purpose : Initialization of the grid coarsening module 
39      !!               1. Read namelist
40      !!               X2. MOVE TO crs_dom.F90 Set the domain definitions for coarse grid
41      !!                 a. Define the coarse grid starting/ending indices on parent grid
42      !!                    Here is where the T-pivot or F-pivot grids are discerned
43      !!                 b. TODO.  Include option for north-centric or equator-centric binning.
44      !!                 (centered only for odd reduction factors; even reduction bins bias equator to the south)
45      !!               3. Mask and mesh creation. => calls to crsfun
46      !!                  a. Use crsfun_mask to generate tmask,umask, vmask, fmask.
47      !!                  b. Use crsfun_coordinates to get coordinates
48      !!                  c. Use crsfun_UV to get horizontal scale factors
49      !!                  d. Use crsfun_TW to get initial vertical scale factors   
50      !!               4. Volumes and weights jes.... TODO. Updates for vvl? Where to do this? crsstp.F90?
51      !!                  a. Calculate initial coarse grid box volumes: pvol_T, pvol_W
52      !!                  b. Calculate initial coarse grid surface-averaging weights
53      !!                  c. Calculate initial coarse grid volume-averaging weights
54      !!                 
55      !!               X5. MOVE TO crs_dom_wri.F90 Using iom_rstput output the initial meshmask.
56      !!               ?. Another set of "masks" to generate
57      !!                  are the u- and v- surface areas for U- and V- area-weighted means.
58      !!                  Need to put this somewhere in section 3?
59      !! jes. What do to about the vvl?  GM.  could separate the weighting (denominator), so
60      !! output C*dA or C*dV as summation not mran, then do mean (division) at moment of output.
61      !! As is, crsfun takes into account vvl.   
62      !!      Talked about pre-setting the surface array to avoid IF/ENDIFS and division.
63      !!      But have then to make that preset array here and elsewhere.
64      !!      that is called every timestep...
65      !!
66      !!               - Read in pertinent data ?
67      !!-------------------------------------------------------------------
68      !! Local variables
[7320]69      INTEGER  :: ji,jj,jk           ! dummy indices
70      INTEGER  :: ierr               ! allocation error status
71      INTEGER  :: ios                ! Local integer output status for namelist read
72      INTEGER  :: iif, iil, ijf, ijl ! indices for tmask_i_crs construction
[5105]73      REAL(wp) :: zmin,zmax
[4015]74      REAL(wp), DIMENSION(:,:,:), POINTER :: zfse3t, zfse3u, zfse3v, zfse3w
75
[7217]76      NAMELIST/namcrs/ nn_factx, nn_facty, nn_msh_crs, nn_crs_kz, ln_crs_wn, ln_crs_top
[4015]77      !!----------------------------------------------------------------------
78      !
79      IF( nn_timing == 1 )  CALL timing_start('crs_init')
80      !
81      IF(lwp) THEN
82         WRITE(numout,*)
83         WRITE(numout,*) 'crs_init : Initializing the grid coarsening module '
84      ENDIF
85
86     !---------------------------------------------------------
87     ! 1. Read Namelist file
88     !---------------------------------------------------------
89     !
[4154]90
[7256]91     REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
92     READ  ( numnam_ref, namcrs, IOSTAT = ios, ERR = 901)
93901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcrs in reference namelist', lwp )
[4154]94
[7256]95     REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
96     READ  ( numnam_cfg, namcrs, IOSTAT = ios, ERR = 902 )
97902  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcrs in configuration namelist', lwp )
98     IF(lwm) WRITE ( numond, namcrs )
[4154]99
[7256]100     IF( .NOT. lk_crs )ln_crs_top = .FALSE. 
101 
[4015]102     IF(lwp) THEN
103        WRITE(numout,*)
104        WRITE(numout,*) 'crs_init: Namelist namcrs '
105        WRITE(numout,*) '   coarsening factor in i-direction      nn_factx   = ', nn_factx
106        WRITE(numout,*) '   coarsening factor in j-direction      nn_facty   = ', nn_facty
107        WRITE(numout,*) '   create (=1) a mesh file or not (=0)   nn_msh_crs = ', nn_msh_crs
108        WRITE(numout,*) '   type of Kz coarsening (0,1,2)         nn_crs_kz  = ', nn_crs_kz
[6797]109
[7256]110        IF( ln_crs_wn )THEN
111           WRITE(numout,*) '   vertical velocities are coarsened '
112        ELSE
113           WRITE(numout,*) '   computed using hdivn '
114        ENDIF
115
116        IF( .NOT. lk_crs )ln_crs_top = .FALSE. 
117
118        IF( ln_crs_top )THEN ; WRITE(numout,*) '   coarsning of physics activated for outputs and BGC model'
119        ELSE                 ; WRITE(numout,*) '   coarsning of physics activated only for outputs'   
120        ENDIF
121
[6797]122        SELECT CASE ( nn_crs_kz )
123           CASE ( 0 ) ; WRITE(numout,*) '   coarsene KZ with MEAN(KZ)'
124           CASE ( 1 ) ; WRITE(numout,*) '   coarsene KZ with MIN(KZ)'
125           CASE ( 2 ) ; WRITE(numout,*) '   coarsene KZ with MAX(KZ)'
126           CASE ( 3 ) ; WRITE(numout,*) '   coarsene KZ with MEANLOG(KZ)'
127           CASE ( 4 ) ; WRITE(numout,*) '   coarsene KZ with MEDIANE(KZ)'
128       END SELECT
[7256]129
[4015]130     ENDIF
131             
132     rfactx_r = 1. / nn_factx
133     rfacty_r = 1. / nn_facty
134
[6772]135
[4015]136     !---------------------------------------------------------
137     ! 2. Define Global Dimensions of the coarsened grid
138     !---------------------------------------------------------
139     CALL crs_dom_def     
140
141     !---------------------------------------------------------
142     ! 3. Mask and Mesh
143     !---------------------------------------------------------
144
145     !  3.a. Get the masks   
146 
147     CALL crs_dom_msk
148
149     !  3.b. Get the coordinates
150     !      Odd-numbered reduction factor, center coordinate on T-cell
151     !      Even-numbered reduction factor, center coordinate on U-,V- faces or f-corner.
152     !     
153        CALL crs_dom_coordinates( gphit, glamt, 'T', gphit_crs, glamt_crs ) 
154        CALL crs_dom_coordinates( gphiu, glamu, 'U', gphiu_crs, glamu_crs )       
155        CALL crs_dom_coordinates( gphiv, glamv, 'V', gphiv_crs, glamv_crs ) 
156        CALL crs_dom_coordinates( gphif, glamf, 'F', gphif_crs, glamf_crs ) 
157
158     !  3.c. Get the horizontal mesh
159
160     !      3.c.1 Horizontal scale factors
161
162     CALL crs_dom_hgr( e1t, e2t, 'T', e1t_crs, e2t_crs )
163     CALL crs_dom_hgr( e1u, e2u, 'U', e1u_crs, e2u_crs )
164     CALL crs_dom_hgr( e1v, e2v, 'V', e1v_crs, e2v_crs )
165     CALL crs_dom_hgr( e1f, e2f, 'F', e1f_crs, e2f_crs )
166
[5105]167     WHERE(e1t_crs == 0._wp) e1t_crs=r_inf
168     WHERE(e1u_crs == 0._wp) e1u_crs=r_inf
169     WHERE(e1v_crs == 0._wp) e1v_crs=r_inf
170     WHERE(e1f_crs == 0._wp) e1f_crs=r_inf
171     WHERE(e2t_crs == 0._wp) e2t_crs=r_inf
172     WHERE(e2u_crs == 0._wp) e2u_crs=r_inf
173     WHERE(e2v_crs == 0._wp) e2v_crs=r_inf
174     WHERE(e2f_crs == 0._wp) e2f_crs=r_inf
175
[4015]176     e1e2t_crs(:,:) = e1t_crs(:,:) * e2t_crs(:,:)
177     
178     
179     !      3.c.2 Coriolis factor 
180
181      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
182
183      CASE ( 0, 1, 4 )           ! mesh on the sphere
184
185         ff_crs(:,:) = 2. * omega * SIN( rad * gphif_crs(:,:) )
186
187      CASE DEFAULT 
188
189       IF(lwp)    WRITE(numout,*) 'crsini.F90. crs_init. Only jphgr_msh = 0, 1 or 4 supported' 
190 
191      END SELECT 
192
193     !    3.d.1 mbathy ( vertical k-levels of bathymetry )     
[10107]194     CALL crs_dom_bat
[4015]195
[10107]196     !    3.d.1 mikt ( isf, top first level)
197     CALL crs_dom_top
[4015]198     
199     !
200     CALL wrk_alloc(jpi, jpj, jpk, zfse3t, zfse3u, zfse3v, zfse3w )
201     !
[6772]202     zfse3t(:,:,:) = e3t_0(:,:,:) !fse3t(:,:,:)
203     zfse3u(:,:,:) = e3u_0(:,:,:) !fse3u(:,:,:)
204     zfse3v(:,:,:) = e3v_0(:,:,:) !fse3v(:,:,:)
205     zfse3w(:,:,:) = e3w_0(:,:,:) !fse3w(:,:,:)
[4015]206
207     !    3.d.2   Surfaces
[5105]208     e2e3u_crs(:,:,:)=0._wp
209     e2e3u_msk(:,:,:)=0._wp
210     e1e3v_crs(:,:,:)=0._wp
211     e1e3v_msk(:,:,:)=0._wp
[4015]212     CALL crs_dom_sfc( tmask, 'W', e1e2w_crs, e1e2w_msk, p_e1=e1t, p_e2=e2t    )
213     CALL crs_dom_sfc( umask, 'U', e2e3u_crs, e2e3u_msk, p_e2=e2u, p_e3=zfse3u )
214     CALL crs_dom_sfc( vmask, 'V', e1e3v_crs, e1e3v_msk, p_e1=e1v, p_e3=zfse3v )
[6772]215
[5007]216     DO jk=1,jpk
217        DO ji=1,jpi_crs
218           DO jj=1,jpj_crs
219
220              facsurfu(ji,jj,jk) = umask_crs(ji,jj,jk) * e2e3u_msk(ji,jj,jk) 
[5601]221              IF( e2e3u_crs(ji,jj,jk) .NE. 0._wp ) facsurfu(ji,jj,jk) = facsurfu(ji,jj,jk) / e2e3u_crs(ji,jj,jk)
[5007]222
223              facsurfv(ji,jj,jk) = vmask_crs(ji,jj,jk) * e1e3v_msk(ji,jj,jk) 
224              IF( e1e3v_crs(ji,jj,jk) .NE. 0._wp ) facsurfv(ji,jj,jk) = facsurfv(ji,jj,jk) / e1e3v_crs(ji,jj,jk)
225
226           ENDDO
227        ENDDO
228     ENDDO
229
[4015]230     !    3.d.3   Vertical scale factors
231     !
[6772]232     CALL crs_dom_e3( e1t, e2t, zfse3t, p_sfc_3d_crs=e1e2w_crs, cd_type='T', p_mask=tmask, p_e3_crs=e3t_0_crs, p_e3_max_crs=e3t_max_0_crs)
233     CALL crs_dom_e3( e1t, e2t, zfse3w, p_sfc_3d_crs=e1e2w_crs, cd_type='W', p_mask=tmask, p_e3_crs=e3w_0_crs, p_e3_max_crs=e3w_max_0_crs)
234     CALL crs_dom_e3( e1u, e2u, zfse3u, p_sfc_2d_crs=e2u_crs  , cd_type='U', p_mask=umask, p_e3_crs=e3u_0_crs, p_e3_max_crs=e3u_max_0_crs)
235     CALL crs_dom_e3( e1v, e2v, zfse3v, p_sfc_2d_crs=e1v_crs  , cd_type='V', p_mask=vmask, p_e3_crs=e3v_0_crs, p_e3_max_crs=e3v_max_0_crs)
[7207]236
[7215]237     WHERE(e3t_max_0_crs == 0._wp) e3t_max_0_crs=r_inf
238     WHERE(e3u_max_0_crs == 0._wp) e3u_max_0_crs=r_inf
239     WHERE(e3v_max_0_crs == 0._wp) e3v_max_0_crs=r_inf
240     WHERE(e3w_max_0_crs == 0._wp) e3w_max_0_crs=r_inf
[7398]241     DO jk = 1, jpk
242        DO ji = 1, jpi_crs
243           DO jj = 1, jpj_crs
244              IF( e3t_max_0_crs(ji,jj,jk) == 0._wp ) e3t_max_0_crs(ji,jj,jk) = e3t_1d(jk)
245              IF( e3w_max_0_crs(ji,jj,jk) == 0._wp ) e3w_max_0_crs(ji,jj,jk) = e3w_1d(jk)
246              IF( e3u_max_0_crs(ji,jj,jk) == 0._wp ) e3u_max_0_crs(ji,jj,jk) = e3t_1d(jk)
247              IF( e3v_max_0_crs(ji,jj,jk) == 0._wp ) e3v_max_0_crs(ji,jj,jk) = e3t_1d(jk)
248           ENDDO
249        ENDDO
250     ENDDO
[6772]251
252#if defined key_vvl
253     e3t_max_n_crs=e3t_max_0_crs
254     e3u_max_n_crs=e3u_max_0_crs
255     e3v_max_n_crs=e3v_max_0_crs
256     e3w_max_n_crs=e3w_max_0_crs
257#endif
258
259     ht_0_crs(:,:)=0._wp
260     DO jk = 1, jpk
261        ht_0_crs(:,:)=ht_0_crs(:,:)+e3t_0_crs(:,:,jk)*tmask_crs(:,:,jk)
262     ENDDO
263
264#if defined key_vvl
265     e3t_0_crs(:,:,:) = e3t_0_crs(:,:,:) * tmask_crs(:,:,:)
266     e3u_0_crs(:,:,:) = e3u_0_crs(:,:,:) * umask_crs(:,:,:)
267     e3v_0_crs(:,:,:) = e3v_0_crs(:,:,:) * vmask_crs(:,:,:)
268     e3w_0_crs(:,:,:) = e3w_0_crs(:,:,:) * tmask_crs(:,:,:)
269#endif
270
[4015]271     ! Reset 0 to e3t_0 or e3w_0
272     DO jk = 1, jpk
273        DO ji = 1, jpi_crs
274           DO jj = 1, jpj_crs
[6772]275              IF( e3t_0_crs(ji,jj,jk) == 0._wp ) e3t_0_crs(ji,jj,jk) = e3t_1d(jk)
276              IF( e3w_0_crs(ji,jj,jk) == 0._wp ) e3w_0_crs(ji,jj,jk) = e3w_1d(jk)
277              IF( e3u_0_crs(ji,jj,jk) == 0._wp ) e3u_0_crs(ji,jj,jk) = e3t_1d(jk)
278              IF( e3v_0_crs(ji,jj,jk) == 0._wp ) e3v_0_crs(ji,jj,jk) = e3t_1d(jk)
[4015]279           ENDDO
280        ENDDO
281     ENDDO
282
[6772]283#if defined key_vvl
284     e3t_b_crs(:,:,:) = e3t_0_crs(:,:,:)
285     e3u_b_crs(:,:,:) = e3u_0_crs(:,:,:)
286     e3v_b_crs(:,:,:) = e3v_0_crs(:,:,:)
287     e3w_b_crs(:,:,:) = e3w_0_crs(:,:,:)
288
289     e3t_n_crs(:,:,:) = e3t_0_crs(:,:,:)
290     e3u_n_crs(:,:,:) = e3u_0_crs(:,:,:)
291     e3v_n_crs(:,:,:) = e3v_0_crs(:,:,:)
292     e3w_n_crs(:,:,:) = e3w_0_crs(:,:,:)
293
294     e3t_a_crs(:,:,:) = e3t_0_crs(:,:,:)
295     e3u_a_crs(:,:,:) = e3u_0_crs(:,:,:)
296     e3v_a_crs(:,:,:) = e3v_0_crs(:,:,:)
297     e3w_a_crs(:,:,:) = e3w_0_crs(:,:,:)
298#endif
299
[4015]300     !    3.d.3   Vertical depth (meters)
[5601]301     !cbr: il semblerait que p_e3=... ne soit pas utile ici !!!!!!!!!
[6772]302     CALL crs_dom_ope( gdept_0, 'MAX', 'T', tmask, gdept_0_crs, p_e3=zfse3t, psgn=1.0 ) 
303     CALL crs_dom_ope( gdepw_0, 'MAX', 'W', tmask, gdepw_0_crs, p_e3=zfse3w, psgn=1.0 )
[7398]304
305     DO jk = 1, jpk
306        DO ji = 1, jpi_crs
307           DO jj = 1, jpj_crs
308              IF( gdept_0_crs(ji,jj,jk) .LE. 0._wp ) gdept_0_crs(ji,jj,jk) = gdept_1d(jk)
309              IF( gdepw_0_crs(ji,jj,jk) .LE. 0._wp ) gdepw_0_crs(ji,jj,jk) = gdepw_1d(jk)
310           ENDDO
311        ENDDO
312     ENDDO
313
[6772]314#if defined key_vvl
315     gdept_n_crs(:,:,:) = gdept_0_crs(:,:,:)
316     gdepw_n_crs(:,:,:) = gdepw_0_crs(:,:,:)
317#endif
[4015]318
[6772]319
[4015]320     !---------------------------------------------------------
321     ! 4.  Coarse grid ocean volume and averaging weights
322     !---------------------------------------------------------
323     CALL crs_dom_facvol( tmask, 'T', e1t, e2t, zfse3t, ocean_volume_crs_t, facvol_t )
324     !
325     bt_crs(:,:,:) = ocean_volume_crs_t(:,:,:) * facvol_t(:,:,:)
326     !
327     r1_bt_crs(:,:,:) = 0._wp 
328     WHERE( bt_crs /= 0._wp ) r1_bt_crs(:,:,:) = 1._wp / bt_crs(:,:,:)
329
330     CALL crs_dom_facvol( tmask, 'W', e1t, e2t, zfse3w, ocean_volume_crs_w, facvol_w )
[5601]331
332
333     !---------------------------------------------------------
[7320]334     ! 5.  Build tmask_i_crs ( for crsdomwri and lib_fortran_crs (glob_sum)
[5601]335     !---------------------------------------------------------
[7320]336     tmask_i_crs(:,:) = tmask_crs(:,:,1)
337     iif = jpreci
338     iil = nlci_crs - jpreci + 1
339     ijf = jpreci
340     ijl = nlcj_crs - jprecj + 1
[5601]341
[7320]342     tmask_i_crs( 1:iif ,    :  ) = 0._wp
343     tmask_i_crs(iil:jpi_crs,    :  ) = 0._wp
344     tmask_i_crs(   :   , 1:ijf ) = 0._wp
345     tmask_i_crs(   :   ,ijl:jpj_crs) = 0._wp
[5601]346
[7320]347     tpol_crs(1:jpiglo_crs,:) = 1._wp
348     fpol_crs(1:jpiglo_crs,:) = 1._wp
349     IF( jperio == 3 .OR. jperio == 4 ) THEN
350        tpol_crs(jpiglo_crs/2+1:jpiglo_crs,:) = 0._wp
351        fpol_crs(       1      :jpiglo_crs,:) = 0._wp
352        IF( mjg_crs(nlej_crs) == jpiglo_crs ) THEN
353           DO ji = iif+1, iil-1
354              tmask_i_crs(ji,nlej_crs-1) = tmask_i_crs(ji,nlej_crs-1) &
355              & * tpol_crs(mig_crs(ji),1)
356           ENDDO
357        ENDIF
358     ENDIF
359     IF( jperio == 5 .OR. jperio == 6 ) THEN
360        tpol_crs(      1       :jpiglo_crs,:)=0._wp
361        fpol_crs(jpiglo_crs/2+1:jpiglo_crs,:)=0._wp
362     ENDIF
363
[4015]364     !
365     !---------------------------------------------------------
[5601]366     ! 6.  Write out coarse meshmask  (see OPA_SRC/DOM/domwri.F90 for ideas later)
[4015]367     !---------------------------------------------------------
368
369     IF( nn_msh_crs > 0 ) THEN
370        CALL dom_grid_crs   ! Save the parent grid information  & Switch to coarse grid domain
371        CALL crs_dom_wri     
372        CALL dom_grid_glo   ! Return to parent grid domain
373     ENDIF
[5105]374   
375
376      rhop_crs(:,:,:)=0._wp ; rhd_crs(:,:,:)=0._wp ; rb2_crs(:,:,:)=0._wp
377
378 
[4015]379     !---------------------------------------------------------
380     ! 7. Finish and clean-up
381     !---------------------------------------------------------
382     CALL wrk_dealloc(jpi, jpj, jpk, zfse3t, zfse3u, zfse3v, zfse3w )
383
[5601]384      IF( nn_timing == 1 )  CALL timing_stop('crs_init')
[4015]385
386   END SUBROUTINE crs_init
387   
388   !!======================================================================
389
390END MODULE crsini
Note: See TracBrowser for help on using the repository browser.