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

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

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