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 @ 7256

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

phaze NEMO routines in CRS branch with nemo_v3_6_STABLE branch at rev 7213 (09-09-2016) (merge -r 5519:7213 )

  • Property svn:keywords set to Id
File size: 13.9 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      REAL(wp) :: zmin,zmax
73      REAL(wp), DIMENSION(:,:,:), POINTER :: zfse3t, zfse3u, zfse3v, zfse3w
74
75      NAMELIST/namcrs/ nn_factx, nn_facty, nn_msh_crs, nn_crs_kz, ln_crs_wn, ln_crs_top
76      !!----------------------------------------------------------------------
77      !
78      IF( nn_timing == 1 )  CALL timing_start('crs_init')
79      !
80      IF(lwp) THEN
81         WRITE(numout,*)
82         WRITE(numout,*) 'crs_init : Initializing the grid coarsening module '
83      ENDIF
84
85     !---------------------------------------------------------
86     ! 1. Read Namelist file
87     !---------------------------------------------------------
88     !
89
90     REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
91     READ  ( numnam_ref, namcrs, IOSTAT = ios, ERR = 901)
92901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcrs in reference namelist', lwp )
93
94     REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
95     READ  ( numnam_cfg, namcrs, IOSTAT = ios, ERR = 902 )
96902  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcrs in configuration namelist', lwp )
97     IF(lwm) WRITE ( numond, namcrs )
98
99     IF( .NOT. lk_crs )ln_crs_top = .FALSE. 
100 
101     IF(lwp) THEN
102        WRITE(numout,*)
103        WRITE(numout,*) 'crs_init: Namelist namcrs '
104        WRITE(numout,*) '   coarsening factor in i-direction      nn_factx   = ', nn_factx
105        WRITE(numout,*) '   coarsening factor in j-direction      nn_facty   = ', nn_facty
106        WRITE(numout,*) '   create (=1) a mesh file or not (=0)   nn_msh_crs = ', nn_msh_crs
107        WRITE(numout,*) '   type of Kz coarsening (0,1,2)         nn_crs_kz  = ', nn_crs_kz
108
109        IF( ln_crs_wn )THEN
110           WRITE(numout,*) '   vertical velocities are coarsened '
111        ELSE
112           WRITE(numout,*) '   computed using hdivn '
113        ENDIF
114
115        IF( .NOT. lk_crs )ln_crs_top = .FALSE. 
116
117        IF( ln_crs_top )THEN ; WRITE(numout,*) '   coarsning of physics activated for outputs and BGC model'
118        ELSE                 ; WRITE(numout,*) '   coarsning of physics activated only for outputs'   
119        ENDIF
120
121        SELECT CASE ( nn_crs_kz )
122           CASE ( 0 ) ; WRITE(numout,*) '   coarsene KZ with MEAN(KZ)'
123           CASE ( 1 ) ; WRITE(numout,*) '   coarsene KZ with MIN(KZ)'
124           CASE ( 2 ) ; WRITE(numout,*) '   coarsene KZ with MAX(KZ)'
125           CASE ( 3 ) ; WRITE(numout,*) '   coarsene KZ with MEANLOG(KZ)'
126           CASE ( 4 ) ; WRITE(numout,*) '   coarsene KZ with MEDIANE(KZ)'
127       END SELECT
128
129     ENDIF
130             
131     rfactx_r = 1. / nn_factx
132     rfacty_r = 1. / nn_facty
133
134
135     !---------------------------------------------------------
136     ! 2. Define Global Dimensions of the coarsened grid
137     !---------------------------------------------------------
138     CALL crs_dom_def     
139
140     !---------------------------------------------------------
141     ! 3. Mask and Mesh
142     !---------------------------------------------------------
143
144     !  3.a. Get the masks   
145 
146     CALL crs_dom_msk
147
148     !  3.b. Get the coordinates
149     !      Odd-numbered reduction factor, center coordinate on T-cell
150     !      Even-numbered reduction factor, center coordinate on U-,V- faces or f-corner.
151     !     
152        CALL crs_dom_coordinates( gphit, glamt, 'T', gphit_crs, glamt_crs ) 
153        CALL crs_dom_coordinates( gphiu, glamu, 'U', gphiu_crs, glamu_crs )       
154        CALL crs_dom_coordinates( gphiv, glamv, 'V', gphiv_crs, glamv_crs ) 
155        CALL crs_dom_coordinates( gphif, glamf, 'F', gphif_crs, glamf_crs ) 
156
157     !  3.c. Get the horizontal mesh
158
159     !      3.c.1 Horizontal scale factors
160
161     CALL crs_dom_hgr( e1t, e2t, 'T', e1t_crs, e2t_crs )
162     CALL crs_dom_hgr( e1u, e2u, 'U', e1u_crs, e2u_crs )
163     CALL crs_dom_hgr( e1v, e2v, 'V', e1v_crs, e2v_crs )
164     CALL crs_dom_hgr( e1f, e2f, 'F', e1f_crs, e2f_crs )
165
166     WHERE(e1t_crs == 0._wp) e1t_crs=r_inf
167     WHERE(e1u_crs == 0._wp) e1u_crs=r_inf
168     WHERE(e1v_crs == 0._wp) e1v_crs=r_inf
169     WHERE(e1f_crs == 0._wp) e1f_crs=r_inf
170     WHERE(e2t_crs == 0._wp) e2t_crs=r_inf
171     WHERE(e2u_crs == 0._wp) e2u_crs=r_inf
172     WHERE(e2v_crs == 0._wp) e2v_crs=r_inf
173     WHERE(e2f_crs == 0._wp) e2f_crs=r_inf
174
175     e1e2t_crs(:,:) = e1t_crs(:,:) * e2t_crs(:,:)
176     
177     
178     !      3.c.2 Coriolis factor 
179
180      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
181
182      CASE ( 0, 1, 4 )           ! mesh on the sphere
183
184         ff_crs(:,:) = 2. * omega * SIN( rad * gphif_crs(:,:) )
185
186      CASE DEFAULT 
187
188       IF(lwp)    WRITE(numout,*) 'crsini.F90. crs_init. Only jphgr_msh = 0, 1 or 4 supported' 
189 
190      END SELECT 
191
192     !    3.d.1 mbathy ( vertical k-levels of bathymetry )     
193
194     CALL crs_dom_bat
195     
196     !
197     CALL wrk_alloc(jpi, jpj, jpk, zfse3t, zfse3u, zfse3v, zfse3w )
198     !
199     zfse3t(:,:,:) = e3t_0(:,:,:) !fse3t(:,:,:)
200     zfse3u(:,:,:) = e3u_0(:,:,:) !fse3u(:,:,:)
201     zfse3v(:,:,:) = e3v_0(:,:,:) !fse3v(:,:,:)
202     zfse3w(:,:,:) = e3w_0(:,:,:) !fse3w(:,:,:)
203
204     !    3.d.2   Surfaces
205     e2e3u_crs(:,:,:)=0._wp
206     e2e3u_msk(:,:,:)=0._wp
207     e1e3v_crs(:,:,:)=0._wp
208     e1e3v_msk(:,:,:)=0._wp
209     CALL crs_dom_sfc( tmask, 'W', e1e2w_crs, e1e2w_msk, p_e1=e1t, p_e2=e2t    )
210     CALL crs_dom_sfc( umask, 'U', e2e3u_crs, e2e3u_msk, p_e2=e2u, p_e3=zfse3u )
211     CALL crs_dom_sfc( vmask, 'V', e1e3v_crs, e1e3v_msk, p_e1=e1v, p_e3=zfse3v )
212
213     DO jk=1,jpk
214        DO ji=1,jpi_crs
215           DO jj=1,jpj_crs
216
217              facsurfu(ji,jj,jk) = umask_crs(ji,jj,jk) * e2e3u_msk(ji,jj,jk) 
218              IF( e2e3u_crs(ji,jj,jk) .NE. 0._wp ) facsurfu(ji,jj,jk) = facsurfu(ji,jj,jk) / e2e3u_crs(ji,jj,jk)
219
220              facsurfv(ji,jj,jk) = vmask_crs(ji,jj,jk) * e1e3v_msk(ji,jj,jk) 
221              IF( e1e3v_crs(ji,jj,jk) .NE. 0._wp ) facsurfv(ji,jj,jk) = facsurfv(ji,jj,jk) / e1e3v_crs(ji,jj,jk)
222
223           ENDDO
224        ENDDO
225     ENDDO
226
227     !    3.d.3   Vertical scale factors
228     !
229     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)
230     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)
231     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)
232     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)
233
234     WHERE(e3t_max_0_crs == 0._wp) e3t_max_0_crs=r_inf
235     WHERE(e3u_max_0_crs == 0._wp) e3u_max_0_crs=r_inf
236     WHERE(e3v_max_0_crs == 0._wp) e3v_max_0_crs=r_inf
237     WHERE(e3w_max_0_crs == 0._wp) e3w_max_0_crs=r_inf
238
239#if defined key_vvl
240     e3t_max_n_crs=e3t_max_0_crs
241     e3u_max_n_crs=e3u_max_0_crs
242     e3v_max_n_crs=e3v_max_0_crs
243     e3w_max_n_crs=e3w_max_0_crs
244#endif
245
246     ht_0_crs(:,:)=0._wp
247     DO jk = 1, jpk
248        ht_0_crs(:,:)=ht_0_crs(:,:)+e3t_0_crs(:,:,jk)*tmask_crs(:,:,jk)
249     ENDDO
250
251#if defined key_vvl
252     e3t_0_crs(:,:,:) = e3t_0_crs(:,:,:) * tmask_crs(:,:,:)
253     e3u_0_crs(:,:,:) = e3u_0_crs(:,:,:) * umask_crs(:,:,:)
254     e3v_0_crs(:,:,:) = e3v_0_crs(:,:,:) * vmask_crs(:,:,:)
255     e3w_0_crs(:,:,:) = e3w_0_crs(:,:,:) * tmask_crs(:,:,:)
256#endif
257
258     ! Reset 0 to e3t_0 or e3w_0
259     DO jk = 1, jpk
260        DO ji = 1, jpi_crs
261           DO jj = 1, jpj_crs
262              IF( e3t_0_crs(ji,jj,jk) == 0._wp ) e3t_0_crs(ji,jj,jk) = e3t_1d(jk)
263              IF( e3w_0_crs(ji,jj,jk) == 0._wp ) e3w_0_crs(ji,jj,jk) = e3w_1d(jk)
264              IF( e3u_0_crs(ji,jj,jk) == 0._wp ) e3u_0_crs(ji,jj,jk) = e3t_1d(jk)
265              IF( e3v_0_crs(ji,jj,jk) == 0._wp ) e3v_0_crs(ji,jj,jk) = e3t_1d(jk)
266           ENDDO
267        ENDDO
268     ENDDO
269
270#if defined key_vvl
271     e3t_b_crs(:,:,:) = e3t_0_crs(:,:,:)
272     e3u_b_crs(:,:,:) = e3u_0_crs(:,:,:)
273     e3v_b_crs(:,:,:) = e3v_0_crs(:,:,:)
274     e3w_b_crs(:,:,:) = e3w_0_crs(:,:,:)
275
276     e3t_n_crs(:,:,:) = e3t_0_crs(:,:,:)
277     e3u_n_crs(:,:,:) = e3u_0_crs(:,:,:)
278     e3v_n_crs(:,:,:) = e3v_0_crs(:,:,:)
279     e3w_n_crs(:,:,:) = e3w_0_crs(:,:,:)
280
281     e3t_a_crs(:,:,:) = e3t_0_crs(:,:,:)
282     e3u_a_crs(:,:,:) = e3u_0_crs(:,:,:)
283     e3v_a_crs(:,:,:) = e3v_0_crs(:,:,:)
284     e3w_a_crs(:,:,:) = e3w_0_crs(:,:,:)
285#endif
286
287     !    3.d.3   Vertical depth (meters)
288     !cbr: il semblerait que p_e3=... ne soit pas utile ici !!!!!!!!!
289     CALL crs_dom_ope( gdept_0, 'MAX', 'T', tmask, gdept_0_crs, p_e3=zfse3t, psgn=1.0 ) 
290     CALL crs_dom_ope( gdepw_0, 'MAX', 'W', tmask, gdepw_0_crs, p_e3=zfse3w, psgn=1.0 )
291#if defined key_vvl
292     gdept_n_crs(:,:,:) = gdept_0_crs(:,:,:)
293     gdepw_n_crs(:,:,:) = gdepw_0_crs(:,:,:)
294#endif
295
296
297     !---------------------------------------------------------
298     ! 4.  Coarse grid ocean volume and averaging weights
299     !---------------------------------------------------------
300     CALL crs_dom_facvol( tmask, 'T', e1t, e2t, zfse3t, ocean_volume_crs_t, facvol_t )
301     !
302     bt_crs(:,:,:) = ocean_volume_crs_t(:,:,:) * facvol_t(:,:,:)
303     !
304     r1_bt_crs(:,:,:) = 0._wp 
305     WHERE( bt_crs /= 0._wp ) r1_bt_crs(:,:,:) = 1._wp / bt_crs(:,:,:)
306
307     CALL crs_dom_facvol( tmask, 'W', e1t, e2t, zfse3w, ocean_volume_crs_w, facvol_w )
308
309
310     !---------------------------------------------------------
311     ! 5.  Coarse grid ocean volume and averaging weights
312     !---------------------------------------------------------
313     !CALL dom_grid_crs   ! Save the parent grid information  & Switch to coarse grid domain
314     !CALL ldf_tra_crs_init
315     !CALL dom_grid_glo   ! Return to parent grid domain
316
317
318     !
319     !---------------------------------------------------------
320     ! 6.  Write out coarse meshmask  (see OPA_SRC/DOM/domwri.F90 for ideas later)
321     !---------------------------------------------------------
322
323     IF( nn_msh_crs > 0 ) THEN
324        CALL dom_grid_crs   ! Save the parent grid information  & Switch to coarse grid domain
325        CALL crs_dom_wri     
326        CALL dom_grid_glo   ! Return to parent grid domain
327     ENDIF
328   
329
330      rhop_crs(:,:,:)=0._wp ; rhd_crs(:,:,:)=0._wp ; rb2_crs(:,:,:)=0._wp
331
332 
333     !---------------------------------------------------------
334     ! 7. Finish and clean-up
335     !---------------------------------------------------------
336     CALL wrk_dealloc(jpi, jpj, jpk, zfse3t, zfse3u, zfse3v, zfse3w )
337
338      IF( nn_timing == 1 )  CALL timing_stop('crs_init')
339
340   END SUBROUTINE crs_init
341   
342   !!======================================================================
343
344END MODULE crsini
Note: See TracBrowser for help on using the repository browser.