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

Last change on this file since 5601 was 5601, checked in by cbricaud, 9 years ago

commit changes/bugfix/... for crs ; ok with time-splitting/fixed volume

File size: 13.0 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
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC  crs_init
28
29   !! * Substitutions
30#  include "domzgr_substitute.h90"
31
32CONTAINS
33   
34   SUBROUTINE crs_init 
35      !!-------------------------------------------------------------------
36      !!                     *** SUBROUTINE crs_init
37      !!  ** Purpose : Initialization of the grid coarsening module 
38      !!               1. Read namelist
39      !!               X2. MOVE TO crs_dom.F90 Set the domain definitions for coarse grid
40      !!                 a. Define the coarse grid starting/ending indices on parent grid
41      !!                    Here is where the T-pivot or F-pivot grids are discerned
42      !!                 b. TODO.  Include option for north-centric or equator-centric binning.
43      !!                 (centered only for odd reduction factors; even reduction bins bias equator to the south)
44      !!               3. Mask and mesh creation. => calls to crsfun
45      !!                  a. Use crsfun_mask to generate tmask,umask, vmask, fmask.
46      !!                  b. Use crsfun_coordinates to get coordinates
47      !!                  c. Use crsfun_UV to get horizontal scale factors
48      !!                  d. Use crsfun_TW to get initial vertical scale factors   
49      !!               4. Volumes and weights jes.... TODO. Updates for vvl? Where to do this? crsstp.F90?
50      !!                  a. Calculate initial coarse grid box volumes: pvol_T, pvol_W
51      !!                  b. Calculate initial coarse grid surface-averaging weights
52      !!                  c. Calculate initial coarse grid volume-averaging weights
53      !!                 
54      !!               X5. MOVE TO crs_dom_wri.F90 Using iom_rstput output the initial meshmask.
55      !!               ?. Another set of "masks" to generate
56      !!                  are the u- and v- surface areas for U- and V- area-weighted means.
57      !!                  Need to put this somewhere in section 3?
58      !! jes. What do to about the vvl?  GM.  could separate the weighting (denominator), so
59      !! output C*dA or C*dV as summation not mran, then do mean (division) at moment of output.
60      !! As is, crsfun takes into account vvl.   
61      !!      Talked about pre-setting the surface array to avoid IF/ENDIFS and division.
62      !!      But have then to make that preset array here and elsewhere.
63      !!      that is called every timestep...
64      !!
65      !!               - Read in pertinent data ?
66      !!-------------------------------------------------------------------
67      !! Local variables
68      INTEGER  :: ji,jj,jk      ! dummy indices
69      INTEGER  :: ierr                                ! allocation error status
70      INTEGER  ::   ios                 ! Local integer output status for namelist read
71      REAL(wp) :: zmin,zmax
72      REAL(wp), DIMENSION(:,:,:), POINTER :: zfse3t, zfse3u, zfse3v, zfse3w
73
74      NAMELIST/namcrs/ nn_factx, nn_facty, nn_binref, nn_msh_crs, nn_crs_kz, ln_crs_wn, ln_crs_top
75      !!----------------------------------------------------------------------
76      !
77      IF( nn_timing == 1 )  CALL timing_start('crs_init')
78      !
79      IF(lwp) THEN
80         WRITE(numout,*)
81         WRITE(numout,*) 'crs_init : Initializing the grid coarsening module '
82      ENDIF
83
84     !---------------------------------------------------------
85     ! 1. Read Namelist file
86     !---------------------------------------------------------
87     !
88
89      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
90      READ  ( numnam_ref, namcrs, IOSTAT = ios, ERR = 901)
91901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcrs in reference namelist', lwp )
92
93      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
94      READ  ( numnam_cfg, namcrs, IOSTAT = ios, ERR = 902 )
95902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcrs in configuration namelist', lwp )
96      IF(lwm) WRITE ( numond, namcrs )
97
98     IF(lwp) THEN
99        WRITE(numout,*)
100        WRITE(numout,*) 'crs_init: Namelist namcrs '
101        WRITE(numout,*) '   coarsening factor in i-direction      nn_factx   = ', nn_factx
102        WRITE(numout,*) '   coarsening factor in j-direction      nn_facty   = ', nn_facty
103        WRITE(numout,*) '   bin centering preference              nn_binref  = ', nn_binref
104        WRITE(numout,*) '   create (=1) a mesh file or not (=0)   nn_msh_crs = ', nn_msh_crs
105        WRITE(numout,*) '   type of Kz coarsening (0,1,2)         nn_crs_kz  = ', nn_crs_kz
106        WRITE(numout,*) '   wn coarsened or computed using hdivn  ln_crs_wn  = ', ln_crs_wn
107     ENDIF
108             
109     rfactx_r = 1. / nn_factx
110     rfacty_r = 1. / nn_facty
111
112     !---------------------------------------------------------
113     ! 2. Define Global Dimensions of the coarsened grid
114     !---------------------------------------------------------
115     CALL crs_dom_def     
116
117     !---------------------------------------------------------
118     ! 3. Mask and Mesh
119     !---------------------------------------------------------
120
121     !     Set up the masks and meshes     
122
123     !  3.a. Get the masks   
124 
125     CALL crs_dom_msk
126
127
128     !  3.b. Get the coordinates
129     !      Odd-numbered reduction factor, center coordinate on T-cell
130     !      Even-numbered reduction factor, center coordinate on U-,V- faces or f-corner.
131     !     
132     IF ( nresty /= 0 .AND. nrestx /= 0 ) THEN
133        CALL crs_dom_coordinates( gphit, glamt, 'T', gphit_crs, glamt_crs ) 
134        CALL crs_dom_coordinates( gphiu, glamu, 'U', gphiu_crs, glamu_crs )       
135        CALL crs_dom_coordinates( gphiv, glamv, 'V', gphiv_crs, glamv_crs ) 
136        CALL crs_dom_coordinates( gphif, glamf, 'F', gphif_crs, glamf_crs ) 
137     ELSEIF ( nresty /= 0 .AND. nrestx == 0 ) THEN
138        CALL crs_dom_coordinates( gphiu, glamu, 'T', gphit_crs, glamt_crs )
139        CALL crs_dom_coordinates( gphiu, glamu, 'U', gphiu_crs, glamu_crs )
140        CALL crs_dom_coordinates( gphif, glamf, 'V', gphiv_crs, glamv_crs )
141        CALL crs_dom_coordinates( gphif, glamf, 'F', gphif_crs, glamf_crs )
142     ELSEIF ( nresty == 0 .AND. nrestx /= 0 ) THEN
143        CALL crs_dom_coordinates( gphiv, glamv, 'T', gphit_crs, glamt_crs )
144        CALL crs_dom_coordinates( gphif, glamf, 'U', gphiu_crs, glamu_crs )
145        CALL crs_dom_coordinates( gphiv, glamv, 'V', gphiv_crs, glamv_crs )
146        CALL crs_dom_coordinates( gphif, glamf, 'F', gphif_crs, glamf_crs )
147     ELSE
148        CALL crs_dom_coordinates( gphif, glamf, 'T', gphit_crs, glamt_crs )
149        CALL crs_dom_coordinates( gphif, glamf, 'U', gphiu_crs, glamu_crs )
150        CALL crs_dom_coordinates( gphif, glamf, 'V', gphiv_crs, glamv_crs )
151        CALL crs_dom_coordinates( gphif, glamf, 'F', gphif_crs, glamf_crs )
152     ENDIF
153
154
155     !  3.c. Get the horizontal mesh
156
157     !      3.c.1 Horizontal scale factors
158
159     CALL crs_dom_hgr( e1t, e2t, 'T', e1t_crs, e2t_crs )
160     CALL crs_dom_hgr( e1u, e2u, 'U', e1u_crs, e2u_crs )
161     CALL crs_dom_hgr( e1v, e2v, 'V', e1v_crs, e2v_crs )
162     CALL crs_dom_hgr( e1f, e2f, 'F', e1f_crs, e2f_crs )
163
164     WHERE(e1t_crs == 0._wp) e1t_crs=r_inf
165     WHERE(e1u_crs == 0._wp) e1u_crs=r_inf
166     WHERE(e1v_crs == 0._wp) e1v_crs=r_inf
167     WHERE(e1f_crs == 0._wp) e1f_crs=r_inf
168     WHERE(e2t_crs == 0._wp) e2t_crs=r_inf
169     WHERE(e2u_crs == 0._wp) e2u_crs=r_inf
170     WHERE(e2v_crs == 0._wp) e2v_crs=r_inf
171     WHERE(e2f_crs == 0._wp) e2f_crs=r_inf
172
173     e1e2t_crs(:,:) = e1t_crs(:,:) * e2t_crs(:,:)
174     
175     
176     !      3.c.2 Coriolis factor 
177
178      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
179
180      CASE ( 0, 1, 4 )           ! mesh on the sphere
181
182         ff_crs(:,:) = 2. * omega * SIN( rad * gphif_crs(:,:) )
183
184      CASE DEFAULT 
185
186       IF(lwp)    WRITE(numout,*) 'crsini.F90. crs_init. Only jphgr_msh = 0, 1 or 4 supported' 
187 
188      END SELECT 
189
190     !    3.d.1 mbathy ( vertical k-levels of bathymetry )     
191
192     CALL crs_dom_bat
193     
194     !
195     CALL wrk_alloc(jpi, jpj, jpk, zfse3t, zfse3u, zfse3v, zfse3w )
196     !
197     zfse3t(:,:,:) = fse3t(:,:,:)
198     zfse3u(:,:,:) = fse3u(:,:,:)
199     zfse3v(:,:,:) = fse3v(:,:,:)
200     zfse3w(:,:,:) = fse3w(:,:,:)
201
202     !    3.d.2   Surfaces
203     e2e3u_crs(:,:,:)=0._wp
204     e2e3u_msk(:,:,:)=0._wp
205     e1e3v_crs(:,:,:)=0._wp
206     e1e3v_msk(:,:,:)=0._wp
207     CALL crs_dom_sfc( tmask, 'W', e1e2w_crs, e1e2w_msk, p_e1=e1t, p_e2=e2t    )
208     CALL crs_dom_sfc( umask, 'U', e2e3u_crs, e2e3u_msk, p_e2=e2u, p_e3=zfse3u )
209     CALL crs_dom_sfc( vmask, 'V', e1e3v_crs, e1e3v_msk, p_e1=e1v, p_e3=zfse3v )
210   
211     !cbr facsurfu(:,:,:) = umask_crs(:,:,:) * e2e3u_msk(:,:,:) / e2e3u_crs(:,:,:)
212     !cbr facsurfv(:,:,:) = vmask_crs(:,:,:) * e1e3v_msk(:,:,:) / e1e3v_crs(:,:,:)
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   
230 
231     CALL crs_dom_e3( e1t, e2t, zfse3t, e1e2w_crs, 'T', tmask, e3t_crs, e3t_max_crs)
232     CALL crs_dom_e3( e1u, e2u, zfse3u, e2e3u_crs, 'U', umask, e3u_crs, e3u_max_crs)
233     CALL crs_dom_e3( e1v, e2v, zfse3v, e1e3v_crs, 'V', vmask, e3v_crs, e3v_max_crs)
234     CALL crs_dom_e3( e1t, e2t, zfse3w, e1e2w_crs, 'W', tmask, e3w_crs, e3w_max_crs)
235     WHERE(e3t_max_crs == 0._wp) e3t_max_crs=r_inf
236     WHERE(e3u_max_crs == 0._wp) e3u_max_crs=r_inf
237     WHERE(e3v_max_crs == 0._wp) e3v_max_crs=r_inf
238     WHERE(e3w_max_crs == 0._wp) e3w_max_crs=r_inf
239
240     ! Reset 0 to e3t_0 or e3w_0
241     DO jk = 1, jpk
242        DO ji = 1, jpi_crs
243           DO jj = 1, jpj_crs
244              IF( e3t_crs(ji,jj,jk) == 0._wp ) e3t_crs(ji,jj,jk) = e3t_1d(jk)
245              IF( e3w_crs(ji,jj,jk) == 0._wp ) e3w_crs(ji,jj,jk) = e3w_1d(jk)
246              IF( e3u_crs(ji,jj,jk) == 0._wp ) e3u_crs(ji,jj,jk) = e3t_1d(jk)
247              IF( e3v_crs(ji,jj,jk) == 0._wp ) e3v_crs(ji,jj,jk) = e3t_1d(jk)
248           ENDDO
249        ENDDO
250     ENDDO
251
252     !    3.d.3   Vertical depth (meters)
253     !cbr: il semblerait que p_e3=... ne soit pas utile ici !!!!!!!!!
254     CALL crs_dom_ope( gdept_0, 'MAX', 'T', tmask, gdept_crs, p_e3=zfse3t, psgn=1.0 ) 
255     CALL crs_dom_ope( gdepw_0, 'MAX', 'W', tmask, gdepw_crs, p_e3=zfse3w, psgn=1.0 )
256
257
258     !---------------------------------------------------------
259     ! 4.  Coarse grid ocean volume and averaging weights
260     !---------------------------------------------------------
261     CALL crs_dom_facvol( tmask, 'T', e1t, e2t, zfse3t, ocean_volume_crs_t, facvol_t )
262     !
263     bt_crs(:,:,:) = ocean_volume_crs_t(:,:,:) * facvol_t(:,:,:)
264     !
265     r1_bt_crs(:,:,:) = 0._wp 
266     WHERE( bt_crs /= 0._wp ) r1_bt_crs(:,:,:) = 1._wp / bt_crs(:,:,:)
267
268     CALL crs_dom_facvol( tmask, 'W', e1t, e2t, zfse3w, ocean_volume_crs_w, facvol_w )
269
270
271     !---------------------------------------------------------
272     ! 5.  Coarse grid ocean volume and averaging weights
273     !---------------------------------------------------------
274     !CALL dom_grid_crs   ! Save the parent grid information  & Switch to coarse grid domain
275     !CALL ldf_tra_crs_init
276     !CALL dom_grid_glo   ! Return to parent grid domain
277
278
279     !
280     !---------------------------------------------------------
281     ! 6.  Write out coarse meshmask  (see OPA_SRC/DOM/domwri.F90 for ideas later)
282     !---------------------------------------------------------
283
284     IF( nn_msh_crs > 0 ) THEN
285        CALL dom_grid_crs   ! Save the parent grid information  & Switch to coarse grid domain
286        CALL crs_dom_wri     
287        CALL dom_grid_glo   ! Return to parent grid domain
288     ENDIF
289   
290
291      rhop_crs(:,:,:)=0._wp ; rhd_crs(:,:,:)=0._wp ; rb2_crs(:,:,:)=0._wp
292
293 
294     !---------------------------------------------------------
295     ! 7. Finish and clean-up
296     !---------------------------------------------------------
297     CALL wrk_dealloc(jpi, jpj, jpk, zfse3t, zfse3u, zfse3v, zfse3w )
298
299      IF( nn_timing == 1 )  CALL timing_stop('crs_init')
300
301   END SUBROUTINE crs_init
302   
303   !!======================================================================
304
305END MODULE crsini
Note: See TracBrowser for help on using the repository browser.