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/UKMO/dev_r10171_test_crs_AMM7/NEMOGCM/NEMO/OPA_SRC/CRS – NEMO

source: branches/UKMO/dev_r10171_test_crs_AMM7/NEMOGCM/NEMO/OPA_SRC/CRS/crsini.F90 @ 10207

Last change on this file since 10207 was 10207, checked in by cmao, 6 years ago

remove svn keyword

File size: 15.6 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_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,*) '   create (=1) a mesh file or not (=0)   nn_msh_crs = ', nn_msh_crs
106        WRITE(numout,*) '   type of Kz coarsening (0,1,2)         nn_crs_kz  = ', nn_crs_kz
107
108        IF( ln_crs_wn )THEN
109           WRITE(numout,*) '   vertical velocities are coarsened '
110        ELSE
111           WRITE(numout,*) '   computed using hdivn '
112        ENDIF
113
114        IF( .NOT. lk_crs )ln_crs_top = .FALSE. 
115
116        IF( ln_crs_top )THEN ; WRITE(numout,*) '   coarsning of physics activated for outputs and BGC model'
117        ELSE                 ; WRITE(numout,*) '   coarsning of physics activated only for outputs'   
118        ENDIF
119
120        SELECT CASE ( nn_crs_kz )
121           CASE ( 0 ) ; WRITE(numout,*) '   coarsene KZ with MEAN(KZ)'
122           CASE ( 1 ) ; WRITE(numout,*) '   coarsene KZ with MIN(KZ)'
123           CASE ( 2 ) ; WRITE(numout,*) '   coarsene KZ with MAX(KZ)'
124           CASE ( 3 ) ; WRITE(numout,*) '   coarsene KZ with MEANLOG(KZ)'
125           CASE ( 4 ) ; WRITE(numout,*) '   coarsene KZ with MEDIANE(KZ)'
126        END SELECT
127
128        WRITE(numout,*) ' '
129        WRITE(numout,*) 'Coarsening factor in i-direction and j-direction: '
130        WRITE(numout,*) 'Be carefull: only nn_factx=nn_facty=3 available'
131     ENDIF
132
133             
134     rfactx_r = 1. / nn_factx
135     rfacty_r = 1. / nn_facty
136
137
138     !---------------------------------------------------------
139     ! 2. Define Global Dimensions of the coarsened grid
140     !---------------------------------------------------------
141     CALL crs_dom_def     
142
143     !---------------------------------------------------------
144     ! 3. Mask and Mesh
145     !---------------------------------------------------------
146
147     !  3.a. Get the masks   
148 
149     CALL crs_dom_msk
150
151     !  3.b. Get the coordinates
152     !      Odd-numbered reduction factor, center coordinate on T-cell
153     !      Even-numbered reduction factor, center coordinate on U-,V- faces or f-corner.
154     !     
155        CALL crs_dom_coordinates( gphit, glamt, 'T', gphit_crs, glamt_crs ) 
156        CALL crs_dom_coordinates( gphiu, glamu, 'U', gphiu_crs, glamu_crs )       
157        CALL crs_dom_coordinates( gphiv, glamv, 'V', gphiv_crs, glamv_crs ) 
158        CALL crs_dom_coordinates( gphif, glamf, 'F', gphif_crs, glamf_crs ) 
159
160     !  3.c. Get the horizontal mesh
161
162     !      3.c.1 Horizontal scale factors
163
164     CALL crs_dom_hgr( e1t, e2t, 'T', e1t_crs, e2t_crs )
165     CALL crs_dom_hgr( e1u, e2u, 'U', e1u_crs, e2u_crs )
166     CALL crs_dom_hgr( e1v, e2v, 'V', e1v_crs, e2v_crs )
167     CALL crs_dom_hgr( e1f, e2f, 'F', e1f_crs, e2f_crs )
168
169     WHERE(e1t_crs == 0._wp) e1t_crs=r_inf
170     WHERE(e1u_crs == 0._wp) e1u_crs=r_inf
171     WHERE(e1v_crs == 0._wp) e1v_crs=r_inf
172     WHERE(e1f_crs == 0._wp) e1f_crs=r_inf
173     WHERE(e2t_crs == 0._wp) e2t_crs=r_inf
174     WHERE(e2u_crs == 0._wp) e2u_crs=r_inf
175     WHERE(e2v_crs == 0._wp) e2v_crs=r_inf
176     WHERE(e2f_crs == 0._wp) e2f_crs=r_inf
177
178     e1e2t_crs(:,:) = e1t_crs(:,:) * e2t_crs(:,:)
179     
180     
181     !      3.c.2 Coriolis factor 
182
183      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
184
185      CASE ( 0, 1, 4 )           ! mesh on the sphere
186
187         ff_crs(:,:) = 2. * omega * SIN( rad * gphif_crs(:,:) )
188
189      CASE DEFAULT 
190
191       IF(lwp)    WRITE(numout,*) 'crsini.F90. crs_init. Only jphgr_msh = 0, 1 or 4 supported' 
192 
193      END SELECT 
194
195     !    3.d.1 mbathy ( vertical k-levels of bathymetry )     
196     CALL crs_dom_bat
197
198     !    3.d.1 mikt ( isf, top first level)
199     CALL crs_dom_top
200     
201     !
202     CALL wrk_alloc(jpi, jpj, jpk, zfse3t, zfse3u, zfse3v, zfse3w )
203     !
204     zfse3t(:,:,:) = e3t_0(:,:,:) !fse3t(:,:,:)
205     zfse3u(:,:,:) = e3u_0(:,:,:) !fse3u(:,:,:)
206     zfse3v(:,:,:) = e3v_0(:,:,:) !fse3v(:,:,:)
207     zfse3w(:,:,:) = e3w_0(:,:,:) !fse3w(:,:,:)
208
209     !    3.d.2   Surfaces
210     e2e3u_crs(:,:,:)=0._wp
211     e2e3u_msk(:,:,:)=0._wp
212     e1e3v_crs(:,:,:)=0._wp
213     e1e3v_msk(:,:,:)=0._wp
214     CALL crs_dom_sfc( tmask, 'W', e1e2w_crs, e1e2w_msk, p_e1=e1t, p_e2=e2t    )
215     CALL crs_dom_sfc( umask, 'U', e2e3u_crs, e2e3u_msk, p_e2=e2u, p_e3=zfse3u )
216     CALL crs_dom_sfc( vmask, 'V', e1e3v_crs, e1e3v_msk, p_e1=e1v, p_e3=zfse3v )
217
218     DO jk=1,jpk
219        DO ji=1,jpi_crs
220           DO jj=1,jpj_crs
221
222              facsurfu(ji,jj,jk) = umask_crs(ji,jj,jk) * e2e3u_msk(ji,jj,jk) 
223              IF( e2e3u_crs(ji,jj,jk) .NE. 0._wp ) facsurfu(ji,jj,jk) = facsurfu(ji,jj,jk) / e2e3u_crs(ji,jj,jk)
224
225              facsurfv(ji,jj,jk) = vmask_crs(ji,jj,jk) * e1e3v_msk(ji,jj,jk) 
226              IF( e1e3v_crs(ji,jj,jk) .NE. 0._wp ) facsurfv(ji,jj,jk) = facsurfv(ji,jj,jk) / e1e3v_crs(ji,jj,jk)
227
228           ENDDO
229        ENDDO
230     ENDDO
231
232     !    3.d.3   Vertical scale factors
233     !
234     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)
235     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)
236     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)
237     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)
238
239     WHERE(e3t_max_0_crs == 0._wp) e3t_max_0_crs=r_inf
240     WHERE(e3u_max_0_crs == 0._wp) e3u_max_0_crs=r_inf
241     WHERE(e3v_max_0_crs == 0._wp) e3v_max_0_crs=r_inf
242     WHERE(e3w_max_0_crs == 0._wp) e3w_max_0_crs=r_inf
243     DO jk = 1, jpk
244        DO ji = 1, jpi_crs
245           DO jj = 1, jpj_crs
246              IF( e3t_max_0_crs(ji,jj,jk) == 0._wp ) e3t_max_0_crs(ji,jj,jk) = e3t_1d(jk)
247              IF( e3w_max_0_crs(ji,jj,jk) == 0._wp ) e3w_max_0_crs(ji,jj,jk) = e3w_1d(jk)
248              IF( e3u_max_0_crs(ji,jj,jk) == 0._wp ) e3u_max_0_crs(ji,jj,jk) = e3t_1d(jk)
249              IF( e3v_max_0_crs(ji,jj,jk) == 0._wp ) e3v_max_0_crs(ji,jj,jk) = e3t_1d(jk)
250           ENDDO
251        ENDDO
252     ENDDO
253
254#if defined key_vvl
255     e3t_max_n_crs=e3t_max_0_crs
256     e3u_max_n_crs=e3u_max_0_crs
257     e3v_max_n_crs=e3v_max_0_crs
258     e3w_max_n_crs=e3w_max_0_crs
259#endif
260
261     ht_0_crs(:,:)=0._wp
262     DO jk = 1, jpk
263        ht_0_crs(:,:)=ht_0_crs(:,:)+e3t_0_crs(:,:,jk)*tmask_crs(:,:,jk)
264     ENDDO
265
266#if defined key_vvl
267     e3t_0_crs(:,:,:) = e3t_0_crs(:,:,:) * tmask_crs(:,:,:)
268     e3u_0_crs(:,:,:) = e3u_0_crs(:,:,:) * umask_crs(:,:,:)
269     e3v_0_crs(:,:,:) = e3v_0_crs(:,:,:) * vmask_crs(:,:,:)
270     e3w_0_crs(:,:,:) = e3w_0_crs(:,:,:) * tmask_crs(:,:,:)
271#endif
272
273     ! Reset 0 to e3t_0 or e3w_0
274     DO jk = 1, jpk
275        DO ji = 1, jpi_crs
276           DO jj = 1, jpj_crs
277              IF( e3t_0_crs(ji,jj,jk) == 0._wp ) e3t_0_crs(ji,jj,jk) = e3t_1d(jk)
278              IF( e3w_0_crs(ji,jj,jk) == 0._wp ) e3w_0_crs(ji,jj,jk) = e3w_1d(jk)
279              IF( e3u_0_crs(ji,jj,jk) == 0._wp ) e3u_0_crs(ji,jj,jk) = e3t_1d(jk)
280              IF( e3v_0_crs(ji,jj,jk) == 0._wp ) e3v_0_crs(ji,jj,jk) = e3t_1d(jk)
281           ENDDO
282        ENDDO
283     ENDDO
284
285#if defined key_vvl
286     e3t_b_crs(:,:,:) = e3t_0_crs(:,:,:)
287     e3u_b_crs(:,:,:) = e3u_0_crs(:,:,:)
288     e3v_b_crs(:,:,:) = e3v_0_crs(:,:,:)
289     e3w_b_crs(:,:,:) = e3w_0_crs(:,:,:)
290
291     e3t_n_crs(:,:,:) = e3t_0_crs(:,:,:)
292     e3u_n_crs(:,:,:) = e3u_0_crs(:,:,:)
293     e3v_n_crs(:,:,:) = e3v_0_crs(:,:,:)
294     e3w_n_crs(:,:,:) = e3w_0_crs(:,:,:)
295
296     e3t_a_crs(:,:,:) = e3t_0_crs(:,:,:)
297     e3u_a_crs(:,:,:) = e3u_0_crs(:,:,:)
298     e3v_a_crs(:,:,:) = e3v_0_crs(:,:,:)
299     e3w_a_crs(:,:,:) = e3w_0_crs(:,:,:)
300#endif
301
302     !    3.d.3   Vertical depth (meters)
303     !cbr: il semblerait que p_e3=... ne soit pas utile ici !!!!!!!!!
304     CALL crs_dom_ope( gdept_0, 'MAX', 'T', tmask, gdept_0_crs, p_e3=zfse3t, psgn=1.0 ) 
305     CALL crs_dom_ope( gdepw_0, 'MAX', 'W', tmask, gdepw_0_crs, p_e3=zfse3w, psgn=1.0 )
306
307     DO jk = 1, jpk
308        DO ji = 1, jpi_crs
309           DO jj = 1, jpj_crs
310              IF( gdept_0_crs(ji,jj,jk) .LE. 0._wp ) gdept_0_crs(ji,jj,jk) = gdept_1d(jk)
311              IF( gdepw_0_crs(ji,jj,jk) .LE. 0._wp ) gdepw_0_crs(ji,jj,jk) = gdepw_1d(jk)
312           ENDDO
313        ENDDO
314     ENDDO
315
316#if defined key_vvl
317     gdept_n_crs(:,:,:) = gdept_0_crs(:,:,:)
318     gdepw_n_crs(:,:,:) = gdepw_0_crs(:,:,:)
319#endif
320
321
322     !---------------------------------------------------------
323     ! 4.  Coarse grid ocean volume and averaging weights
324     !---------------------------------------------------------
325     CALL crs_dom_facvol( tmask, 'T', e1t, e2t, zfse3t, ocean_volume_crs_t, facvol_t )
326     !
327     bt_crs(:,:,:) = ocean_volume_crs_t(:,:,:) * facvol_t(:,:,:)
328     !
329     r1_bt_crs(:,:,:) = 0._wp 
330     WHERE( bt_crs /= 0._wp ) r1_bt_crs(:,:,:) = 1._wp / bt_crs(:,:,:)
331
332     CALL crs_dom_facvol( tmask, 'W', e1t, e2t, zfse3w, ocean_volume_crs_w, facvol_w )
333
334
335     !---------------------------------------------------------
336     ! 5.  Build tmask_i_crs ( for crsdomwri and lib_fortran_crs (glob_sum)
337     !---------------------------------------------------------
338     tmask_i_crs(:,:) = tmask_crs(:,:,1)
339     iif = jpreci
340     iil = nlci_crs - jpreci + 1
341     ijf = jpreci
342     ijl = nlcj_crs - jprecj + 1
343
344     tmask_i_crs( 1:iif ,    :  ) = 0._wp
345     tmask_i_crs(iil:jpi_crs,    :  ) = 0._wp
346     tmask_i_crs(   :   , 1:ijf ) = 0._wp
347     tmask_i_crs(   :   ,ijl:jpj_crs) = 0._wp
348
349     tpol_crs(1:jpiglo_crs,:) = 1._wp
350     fpol_crs(1:jpiglo_crs,:) = 1._wp
351     IF( jperio == 3 .OR. jperio == 4 ) THEN
352        tpol_crs(jpiglo_crs/2+1:jpiglo_crs,:) = 0._wp
353        fpol_crs(       1      :jpiglo_crs,:) = 0._wp
354        IF( mjg_crs(nlej_crs) == jpiglo_crs ) THEN
355           DO ji = iif+1, iil-1
356              tmask_i_crs(ji,nlej_crs-1) = tmask_i_crs(ji,nlej_crs-1) &
357              & * tpol_crs(mig_crs(ji),1)
358           ENDDO
359        ENDIF
360     ENDIF
361     IF( jperio == 5 .OR. jperio == 6 ) THEN
362        tpol_crs(      1       :jpiglo_crs,:)=0._wp
363        fpol_crs(jpiglo_crs/2+1:jpiglo_crs,:)=0._wp
364     ENDIF
365
366     !
367     !---------------------------------------------------------
368     ! 6.  Write out coarse meshmask  (see OPA_SRC/DOM/domwri.F90 for ideas later)
369     !---------------------------------------------------------
370
371     IF( nn_msh_crs > 0 ) THEN
372        CALL dom_grid_crs   ! Save the parent grid information  & Switch to coarse grid domain
373        CALL crs_dom_wri     
374        CALL dom_grid_glo   ! Return to parent grid domain
375     ENDIF
376   
377
378      rhop_crs(:,:,:)=0._wp ; rhd_crs(:,:,:)=0._wp ; rb2_crs(:,:,:)=0._wp
379
380 
381     !---------------------------------------------------------
382     ! 7. Finish and clean-up
383     !---------------------------------------------------------
384     CALL wrk_dealloc(jpi, jpj, jpk, zfse3t, zfse3u, zfse3v, zfse3w )
385
386      IF( nn_timing == 1 )  CALL timing_stop('crs_init')
387
388   END SUBROUTINE crs_init
389   
390   !!======================================================================
391
392END MODULE crsini
Note: See TracBrowser for help on using the repository browser.