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

source: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsini.F90 @ 7795

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

coarsening branch: first implementation of coarsening in PISCES
  • Property svn:keywords set to Id
File size: 15.5 KB
Line 
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
22   USE ldftra_crs
23   USE ieee_arithmetic
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
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
73      REAL(wp) :: zmin,zmax
74      REAL(wp), DIMENSION(:,:,:), POINTER :: zfse3t, zfse3u, zfse3v, zfse3w
75
76      NAMELIST/namcrs/ nn_factx, nn_facty, nn_msh_crs, nn_crs_kz, ln_crs_wn, ln_crs_top
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     !
90
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 )
94
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 )
99
100     IF( .NOT. lk_crs )ln_crs_top = .FALSE. 
101 
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
109
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
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
129
130     ENDIF
131             
132     rfactx_r = 1. / nn_factx
133     rfacty_r = 1. / nn_facty
134
135
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
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
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 )     
194
195     CALL crs_dom_bat
196     
197     !
198     CALL wrk_alloc(jpi, jpj, jpk, zfse3t, zfse3u, zfse3v, zfse3w )
199     !
200     zfse3t(:,:,:) = e3t_0(:,:,:) !fse3t(:,:,:)
201     zfse3u(:,:,:) = e3u_0(:,:,:) !fse3u(:,:,:)
202     zfse3v(:,:,:) = e3v_0(:,:,:) !fse3v(:,:,:)
203     zfse3w(:,:,:) = e3w_0(:,:,:) !fse3w(:,:,:)
204
205     !    3.d.2   Surfaces
206     e2e3u_crs(:,:,:)=0._wp
207     e2e3u_msk(:,:,:)=0._wp
208     e1e3v_crs(:,:,:)=0._wp
209     e1e3v_msk(:,:,:)=0._wp
210     CALL crs_dom_sfc( tmask, 'W', e1e2w_crs, e1e2w_msk, p_e1=e1t, p_e2=e2t    )
211     CALL crs_dom_sfc( umask, 'U', e2e3u_crs, e2e3u_msk, p_e2=e2u, p_e3=zfse3u )
212     CALL crs_dom_sfc( vmask, 'V', e1e3v_crs, e1e3v_msk, p_e1=e1v, p_e3=zfse3v )
213
214     DO jk=1,jpk
215        DO ji=1,jpi_crs
216           DO jj=1,jpj_crs
217
218              facsurfu(ji,jj,jk) = umask_crs(ji,jj,jk) * e2e3u_msk(ji,jj,jk) 
219              IF( e2e3u_crs(ji,jj,jk) .NE. 0._wp ) facsurfu(ji,jj,jk) = facsurfu(ji,jj,jk) / e2e3u_crs(ji,jj,jk)
220
221              facsurfv(ji,jj,jk) = vmask_crs(ji,jj,jk) * e1e3v_msk(ji,jj,jk) 
222              IF( e1e3v_crs(ji,jj,jk) .NE. 0._wp ) facsurfv(ji,jj,jk) = facsurfv(ji,jj,jk) / e1e3v_crs(ji,jj,jk)
223
224           ENDDO
225        ENDDO
226     ENDDO
227
228     !    3.d.3   Vertical scale factors
229     !
230     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)
231     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)
232     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)
233     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)
234
235     WHERE(e3t_max_0_crs == 0._wp) e3t_max_0_crs=r_inf
236     WHERE(e3u_max_0_crs == 0._wp) e3u_max_0_crs=r_inf
237     WHERE(e3v_max_0_crs == 0._wp) e3v_max_0_crs=r_inf
238     WHERE(e3w_max_0_crs == 0._wp) e3w_max_0_crs=r_inf
239     DO jk = 1, jpk
240        DO ji = 1, jpi_crs
241           DO jj = 1, jpj_crs
242              IF( e3t_max_0_crs(ji,jj,jk) == 0._wp ) e3t_max_0_crs(ji,jj,jk) = e3t_1d(jk)
243              IF( e3w_max_0_crs(ji,jj,jk) == 0._wp ) e3w_max_0_crs(ji,jj,jk) = e3w_1d(jk)
244              IF( e3u_max_0_crs(ji,jj,jk) == 0._wp ) e3u_max_0_crs(ji,jj,jk) = e3t_1d(jk)
245              IF( e3v_max_0_crs(ji,jj,jk) == 0._wp ) e3v_max_0_crs(ji,jj,jk) = e3t_1d(jk)
246           ENDDO
247        ENDDO
248     ENDDO
249
250#if defined key_vvl
251     e3t_max_n_crs=e3t_max_0_crs
252     e3u_max_n_crs=e3u_max_0_crs
253     e3v_max_n_crs=e3v_max_0_crs
254     e3w_max_n_crs=e3w_max_0_crs
255#endif
256
257     ht_0_crs(:,:)=0._wp
258     DO jk = 1, jpk
259        ht_0_crs(:,:)=ht_0_crs(:,:)+e3t_0_crs(:,:,jk)*tmask_crs(:,:,jk)
260     ENDDO
261
262#if defined key_vvl
263     e3t_0_crs(:,:,:) = e3t_0_crs(:,:,:) * tmask_crs(:,:,:)
264     e3u_0_crs(:,:,:) = e3u_0_crs(:,:,:) * umask_crs(:,:,:)
265     e3v_0_crs(:,:,:) = e3v_0_crs(:,:,:) * vmask_crs(:,:,:)
266     e3w_0_crs(:,:,:) = e3w_0_crs(:,:,:) * tmask_crs(:,:,:)
267#endif
268
269     ! Reset 0 to e3t_0 or e3w_0
270     DO jk = 1, jpk
271        DO ji = 1, jpi_crs
272           DO jj = 1, jpj_crs
273              IF( e3t_0_crs(ji,jj,jk) == 0._wp ) e3t_0_crs(ji,jj,jk) = e3t_1d(jk)
274              IF( e3w_0_crs(ji,jj,jk) == 0._wp ) e3w_0_crs(ji,jj,jk) = e3w_1d(jk)
275              IF( e3u_0_crs(ji,jj,jk) == 0._wp ) e3u_0_crs(ji,jj,jk) = e3t_1d(jk)
276              IF( e3v_0_crs(ji,jj,jk) == 0._wp ) e3v_0_crs(ji,jj,jk) = e3t_1d(jk)
277           ENDDO
278        ENDDO
279     ENDDO
280
281#if defined key_vvl
282     e3t_b_crs(:,:,:) = e3t_0_crs(:,:,:)
283     e3u_b_crs(:,:,:) = e3u_0_crs(:,:,:)
284     e3v_b_crs(:,:,:) = e3v_0_crs(:,:,:)
285     e3w_b_crs(:,:,:) = e3w_0_crs(:,:,:)
286
287     e3t_n_crs(:,:,:) = e3t_0_crs(:,:,:)
288     e3u_n_crs(:,:,:) = e3u_0_crs(:,:,:)
289     e3v_n_crs(:,:,:) = e3v_0_crs(:,:,:)
290     e3w_n_crs(:,:,:) = e3w_0_crs(:,:,:)
291
292     e3t_a_crs(:,:,:) = e3t_0_crs(:,:,:)
293     e3u_a_crs(:,:,:) = e3u_0_crs(:,:,:)
294     e3v_a_crs(:,:,:) = e3v_0_crs(:,:,:)
295     e3w_a_crs(:,:,:) = e3w_0_crs(:,:,:)
296#endif
297
298     !    3.d.3   Vertical depth (meters)
299     !cbr: il semblerait que p_e3=... ne soit pas utile ici !!!!!!!!!
300     CALL crs_dom_ope( gdept_0, 'MAX', 'T', tmask, gdept_0_crs, p_e3=zfse3t, psgn=1.0 ) 
301     CALL crs_dom_ope( gdepw_0, 'MAX', 'W', tmask, gdepw_0_crs, p_e3=zfse3w, psgn=1.0 )
302
303     DO jk = 1, jpk
304        DO ji = 1, jpi_crs
305           DO jj = 1, jpj_crs
306              IF( gdept_0_crs(ji,jj,jk) .LE. 0._wp ) gdept_0_crs(ji,jj,jk) = gdept_1d(jk)
307              IF( gdepw_0_crs(ji,jj,jk) .LE. 0._wp ) gdepw_0_crs(ji,jj,jk) = gdepw_1d(jk)
308           ENDDO
309        ENDDO
310     ENDDO
311
312#if defined key_vvl
313     gdept_n_crs(:,:,:) = gdept_0_crs(:,:,:)
314     gdepw_n_crs(:,:,:) = gdepw_0_crs(:,:,:)
315#endif
316
317
318     !---------------------------------------------------------
319     ! 4.  Coarse grid ocean volume and averaging weights
320     !---------------------------------------------------------
321     CALL crs_dom_facvol( tmask, 'T', e1t, e2t, zfse3t, ocean_volume_crs_t, facvol_t )
322     !
323     bt_crs(:,:,:) = ocean_volume_crs_t(:,:,:) * facvol_t(:,:,:)
324     !
325     r1_bt_crs(:,:,:) = 0._wp 
326     WHERE( bt_crs /= 0._wp ) r1_bt_crs(:,:,:) = 1._wp / bt_crs(:,:,:)
327
328     CALL crs_dom_facvol( tmask, 'W', e1t, e2t, zfse3w, ocean_volume_crs_w, facvol_w )
329
330
331     !---------------------------------------------------------
332     ! 5.  Build tmask_i_crs ( for crsdomwri and lib_fortran_crs (glob_sum)
333     !---------------------------------------------------------
334     tmask_i_crs(:,:) = tmask_crs(:,:,1)
335     iif = jpreci
336     iil = nlci_crs - jpreci + 1
337     ijf = jpreci
338     ijl = nlcj_crs - jprecj + 1
339
340     tmask_i_crs( 1:iif ,    :  ) = 0._wp
341     tmask_i_crs(iil:jpi_crs,    :  ) = 0._wp
342     tmask_i_crs(   :   , 1:ijf ) = 0._wp
343     tmask_i_crs(   :   ,ijl:jpj_crs) = 0._wp
344
345     tpol_crs(1:jpiglo_crs,:) = 1._wp
346     fpol_crs(1:jpiglo_crs,:) = 1._wp
347     IF( jperio == 3 .OR. jperio == 4 ) THEN
348        tpol_crs(jpiglo_crs/2+1:jpiglo_crs,:) = 0._wp
349        fpol_crs(       1      :jpiglo_crs,:) = 0._wp
350        IF( mjg_crs(nlej_crs) == jpiglo_crs ) THEN
351           DO ji = iif+1, iil-1
352              tmask_i_crs(ji,nlej_crs-1) = tmask_i_crs(ji,nlej_crs-1) &
353              & * tpol_crs(mig_crs(ji),1)
354           ENDDO
355        ENDIF
356     ENDIF
357     IF( jperio == 5 .OR. jperio == 6 ) THEN
358        tpol_crs(      1       :jpiglo_crs,:)=0._wp
359        fpol_crs(jpiglo_crs/2+1:jpiglo_crs,:)=0._wp
360     ENDIF
361
362     !
363     !---------------------------------------------------------
364     ! 6.  Write out coarse meshmask  (see OPA_SRC/DOM/domwri.F90 for ideas later)
365     !---------------------------------------------------------
366
367     IF( nn_msh_crs > 0 ) THEN
368        CALL dom_grid_crs   ! Save the parent grid information  & Switch to coarse grid domain
369        CALL crs_dom_wri     
370        CALL dom_grid_glo   ! Return to parent grid domain
371     ENDIF
372   
373
374      rhop_crs(:,:,:)=0._wp ; rhd_crs(:,:,:)=0._wp ; rb2_crs(:,:,:)=0._wp
375
376 
377     !---------------------------------------------------------
378     ! 7. Finish and clean-up
379     !---------------------------------------------------------
380     CALL wrk_dealloc(jpi, jpj, jpk, zfse3t, zfse3u, zfse3v, zfse3w )
381
382      IF( nn_timing == 1 )  CALL timing_stop('crs_init')
383
384   END SUBROUTINE crs_init
385   
386   !!======================================================================
387
388END MODULE crsini
Note: See TracBrowser for help on using the repository browser.