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

Last change on this file since 6101 was 6101, checked in by cbricaud, 8 years ago

correction of bugs from last update and improvments for CRS
  • Property svn:keywords set to Id
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
32   !! $Id$
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     ENDIF
109             
110     rfactx_r = 1. / nn_factx
111     rfacty_r = 1. / nn_facty
112
113     !---------------------------------------------------------
114     ! 2. Define Global Dimensions of the coarsened grid
115     !---------------------------------------------------------
116     CALL crs_dom_def     
117
118     !---------------------------------------------------------
119     ! 3. Mask and Mesh
120     !---------------------------------------------------------
121
122     !     Set up the masks and meshes     
123
124     !  3.a. Get the masks   
125 
126     CALL crs_dom_msk
127
128
129     !  3.b. Get the coordinates
130     !      Odd-numbered reduction factor, center coordinate on T-cell
131     !      Even-numbered reduction factor, center coordinate on U-,V- faces or f-corner.
132     !     
133     IF ( nresty /= 0 .AND. nrestx /= 0 ) THEN
134        CALL crs_dom_coordinates( gphit, glamt, 'T', gphit_crs, glamt_crs ) 
135        CALL crs_dom_coordinates( gphiu, glamu, 'U', gphiu_crs, glamu_crs )       
136        CALL crs_dom_coordinates( gphiv, glamv, 'V', gphiv_crs, glamv_crs ) 
137        CALL crs_dom_coordinates( gphif, glamf, 'F', gphif_crs, glamf_crs ) 
138     ELSEIF ( nresty /= 0 .AND. nrestx == 0 ) THEN
139        CALL crs_dom_coordinates( gphiu, glamu, 'T', gphit_crs, glamt_crs )
140        CALL crs_dom_coordinates( gphiu, glamu, 'U', gphiu_crs, glamu_crs )
141        CALL crs_dom_coordinates( gphif, glamf, 'V', gphiv_crs, glamv_crs )
142        CALL crs_dom_coordinates( gphif, glamf, 'F', gphif_crs, glamf_crs )
143     ELSEIF ( nresty == 0 .AND. nrestx /= 0 ) THEN
144        CALL crs_dom_coordinates( gphiv, glamv, 'T', gphit_crs, glamt_crs )
145        CALL crs_dom_coordinates( gphif, glamf, 'U', gphiu_crs, glamu_crs )
146        CALL crs_dom_coordinates( gphiv, glamv, 'V', gphiv_crs, glamv_crs )
147        CALL crs_dom_coordinates( gphif, glamf, 'F', gphif_crs, glamf_crs )
148     ELSE
149        CALL crs_dom_coordinates( gphif, glamf, 'T', gphit_crs, glamt_crs )
150        CALL crs_dom_coordinates( gphif, glamf, 'U', gphiu_crs, glamu_crs )
151        CALL crs_dom_coordinates( gphif, glamf, 'V', gphiv_crs, glamv_crs )
152        CALL crs_dom_coordinates( gphif, glamf, 'F', gphif_crs, glamf_crs )
153     ENDIF
154
155
156     !  3.c. Get the horizontal mesh
157
158     !      3.c.1 Horizontal scale factors
159
160     CALL crs_dom_hgr( e1t, e2t, 'T', e1t_crs, e2t_crs )
161     CALL crs_dom_hgr( e1u, e2u, 'U', e1u_crs, e2u_crs )
162     CALL crs_dom_hgr( e1v, e2v, 'V', e1v_crs, e2v_crs )
163     CALL crs_dom_hgr( e1f, e2f, 'F', e1f_crs, e2f_crs )
164
165     WHERE(e1t_crs == 0._wp) e1t_crs=r_inf
166     WHERE(e1u_crs == 0._wp) e1u_crs=r_inf
167     WHERE(e1v_crs == 0._wp) e1v_crs=r_inf
168     WHERE(e1f_crs == 0._wp) e1f_crs=r_inf
169     WHERE(e2t_crs == 0._wp) e2t_crs=r_inf
170     WHERE(e2u_crs == 0._wp) e2u_crs=r_inf
171     WHERE(e2v_crs == 0._wp) e2v_crs=r_inf
172     WHERE(e2f_crs == 0._wp) e2f_crs=r_inf
173
174     e1e2t_crs(:,:) = e1t_crs(:,:) * e2t_crs(:,:)
175     
176     
177     !      3.c.2 Coriolis factor 
178
179      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
180
181      CASE ( 0, 1, 4 )           ! mesh on the sphere
182
183         ff_crs(:,:) = 2. * omega * SIN( rad * gphif_crs(:,:) )
184
185      CASE DEFAULT 
186
187       IF(lwp)    WRITE(numout,*) 'crsini.F90. crs_init. Only jphgr_msh = 0, 1 or 4 supported' 
188 
189      END SELECT 
190
191     !    3.d.1 mbathy ( vertical k-levels of bathymetry )     
192
193     CALL crs_dom_bat
194     
195     !
196     CALL wrk_alloc(jpi, jpj, jpk, zfse3t, zfse3u, zfse3v, zfse3w )
197     !
198     zfse3t(:,:,:) = fse3t(:,:,:)
199     zfse3u(:,:,:) = fse3u(:,:,:)
200     zfse3v(:,:,:) = fse3v(:,:,:)
201     zfse3w(:,:,:) = fse3w(:,:,:)
202
203     !    3.d.2   Surfaces
204     e2e3u_crs(:,:,:)=0._wp
205     e2e3u_msk(:,:,:)=0._wp
206     e1e3v_crs(:,:,:)=0._wp
207     e1e3v_msk(:,:,:)=0._wp
208     CALL crs_dom_sfc( tmask, 'W', e1e2w_crs, e1e2w_msk, p_e1=e1t, p_e2=e2t    )
209     CALL crs_dom_sfc( umask, 'U', e2e3u_crs, e2e3u_msk, p_e2=e2u, p_e3=zfse3u )
210     CALL crs_dom_sfc( vmask, 'V', e1e3v_crs, e1e3v_msk, p_e1=e1v, p_e3=zfse3v )
211   
212     !cbr facsurfu(:,:,:) = umask_crs(:,:,:) * e2e3u_msk(:,:,:) / e2e3u_crs(:,:,:)
213     !cbr facsurfv(:,:,:) = vmask_crs(:,:,:) * e1e3v_msk(:,:,:) / e1e3v_crs(:,:,:)
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   
231 
232     CALL crs_dom_e3( e1t, e2t, zfse3t, e1e2w_crs, 'T', tmask, e3t_crs, e3t_max_crs)
233     CALL crs_dom_e3( e1u, e2u, zfse3u, e2e3u_crs, 'U', umask, e3u_crs, e3u_max_crs)
234     CALL crs_dom_e3( e1v, e2v, zfse3v, e1e3v_crs, 'V', vmask, e3v_crs, e3v_max_crs)
235     CALL crs_dom_e3( e1t, e2t, zfse3w, e1e2w_crs, 'W', tmask, e3w_crs, e3w_max_crs)
236     WHERE(e3t_max_crs == 0._wp) e3t_max_crs=r_inf
237     WHERE(e3u_max_crs == 0._wp) e3u_max_crs=r_inf
238     WHERE(e3v_max_crs == 0._wp) e3v_max_crs=r_inf
239     WHERE(e3w_max_crs == 0._wp) e3w_max_crs=r_inf
240
241     ! Reset 0 to e3t_0 or e3w_0
242     DO jk = 1, jpk
243        DO ji = 1, jpi_crs
244           DO jj = 1, jpj_crs
245              IF( e3t_crs(ji,jj,jk) == 0._wp ) e3t_crs(ji,jj,jk) = e3t_1d(jk)
246              IF( e3w_crs(ji,jj,jk) == 0._wp ) e3w_crs(ji,jj,jk) = e3w_1d(jk)
247              IF( e3u_crs(ji,jj,jk) == 0._wp ) e3u_crs(ji,jj,jk) = e3t_1d(jk)
248              IF( e3v_crs(ji,jj,jk) == 0._wp ) e3v_crs(ji,jj,jk) = e3t_1d(jk)
249           ENDDO
250        ENDDO
251     ENDDO
252
253     !    3.d.3   Vertical depth (meters)
254     !cbr: il semblerait que p_e3=... ne soit pas utile ici !!!!!!!!!
255     CALL crs_dom_ope( gdept_0, 'MAX', 'T', tmask, gdept_crs, p_e3=zfse3t, psgn=1.0 ) 
256     CALL crs_dom_ope( gdepw_0, 'MAX', 'W', tmask, gdepw_crs, p_e3=zfse3w, psgn=1.0 )
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.