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

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

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