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/2013/dev_r3940_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS – NEMO

source: branches/2013/dev_r3940_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsini.F90 @ 4015

Last change on this file since 4015 was 4015, checked in by cetlod, 11 years ago

2013/dev_r3940_CNRS4_IOCRS: 1st step, add new routines for outputs coarsening
File size: 10.8 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      REAL(wp), DIMENSION(:,:,:), POINTER :: zfse3t, zfse3u, zfse3v, zfse3w
70
71      NAMELIST/namcrs/ nn_factx, nn_facty, nn_binref, nn_msh_crs, nn_crs_kz, ln_crs_wn
72
73
74      !!----------------------------------------------------------------------
75      !
76      IF( nn_timing == 1 )  CALL timing_start('crs_init')
77      !
78      IF(lwp) THEN
79         WRITE(numout,*)
80         WRITE(numout,*) 'crs_init : Initializing the grid coarsening module '
81      ENDIF
82
83     !---------------------------------------------------------
84     ! 1. Read Namelist file
85     !---------------------------------------------------------
86     !
87     REWIND( numnam )            ! Read Namelist namtra_qsr : ratio and length of penetration
88     READ  ( numnam, namcrs )        ! Namelist namcrs : Grid-coarsening utility namelist
89     IF(lwp) THEN
90        WRITE(numout,*)
91        WRITE(numout,*) 'crs_init: Namelist namcrs '
92        WRITE(numout,*) '   coarsening factor in i-direction      nn_factx   = ', nn_factx
93        WRITE(numout,*) '   coarsening factor in j-direction      nn_facty   = ', nn_facty
94        WRITE(numout,*) '   bin centering preference              nn_binref  = ', nn_binref
95        WRITE(numout,*) '   create (=1) a mesh file or not (=0)   nn_msh_crs = ', nn_msh_crs
96        WRITE(numout,*) '   type of Kz coarsening (0,1,2)         nn_crs_kz  = ', nn_crs_kz
97        WRITE(numout,*) '   wn coarsened or computed using hdivn  ln_crs_wn  = ', ln_crs_wn
98     ENDIF
99             
100     rfactx_r = 1. / nn_factx
101     rfacty_r = 1. / nn_facty
102
103     !---------------------------------------------------------
104     ! 2. Define Global Dimensions of the coarsened grid
105     !---------------------------------------------------------
106     CALL crs_dom_def     
107
108     !---------------------------------------------------------
109     ! 3. Mask and Mesh
110     !---------------------------------------------------------
111
112     !     Set up the masks and meshes     
113
114     !  3.a. Get the masks   
115 
116     CALL crs_dom_msk
117
118
119     !  3.b. Get the coordinates
120     !      Odd-numbered reduction factor, center coordinate on T-cell
121     !      Even-numbered reduction factor, center coordinate on U-,V- faces or f-corner.
122     !     
123     IF ( nresty /= 0 .AND. nrestx /= 0 ) THEN
124        CALL crs_dom_coordinates( gphit, glamt, 'T', gphit_crs, glamt_crs ) 
125        CALL crs_dom_coordinates( gphiu, glamu, 'U', gphiu_crs, glamu_crs )       
126        CALL crs_dom_coordinates( gphiv, glamv, 'V', gphiv_crs, glamv_crs ) 
127        CALL crs_dom_coordinates( gphif, glamf, 'F', gphif_crs, glamf_crs ) 
128     ELSEIF ( nresty /= 0 .AND. nrestx == 0 ) THEN
129        CALL crs_dom_coordinates( gphiu, glamu, 'T', gphit_crs, glamt_crs )
130        CALL crs_dom_coordinates( gphiu, glamu, 'U', gphiu_crs, glamu_crs )
131        CALL crs_dom_coordinates( gphif, glamf, 'V', gphiv_crs, glamv_crs )
132        CALL crs_dom_coordinates( gphif, glamf, 'F', gphif_crs, glamf_crs )
133     ELSEIF ( nresty == 0 .AND. nrestx /= 0 ) THEN
134        CALL crs_dom_coordinates( gphiv, glamv, 'T', gphit_crs, glamt_crs )
135        CALL crs_dom_coordinates( gphif, glamf, '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     ELSE
139        CALL crs_dom_coordinates( gphif, glamf, 'T', gphit_crs, glamt_crs )
140        CALL crs_dom_coordinates( gphif, glamf, '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     ENDIF
144
145
146     !  3.c. Get the horizontal mesh
147
148     !      3.c.1 Horizontal scale factors
149
150     CALL crs_dom_hgr( e1t, e2t, 'T', e1t_crs, e2t_crs )
151     CALL crs_dom_hgr( e1u, e2u, 'U', e1u_crs, e2u_crs )
152     CALL crs_dom_hgr( e1v, e2v, 'V', e1v_crs, e2v_crs )
153     CALL crs_dom_hgr( e1f, e2f, 'F', e1f_crs, e2f_crs )
154
155     e1e2t_crs(:,:) = e1t_crs(:,:) * e2t_crs(:,:)
156     
157     
158     !      3.c.2 Coriolis factor 
159
160      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
161
162      CASE ( 0, 1, 4 )           ! mesh on the sphere
163
164         ff_crs(:,:) = 2. * omega * SIN( rad * gphif_crs(:,:) )
165
166      CASE DEFAULT 
167
168       IF(lwp)    WRITE(numout,*) 'crsini.F90. crs_init. Only jphgr_msh = 0, 1 or 4 supported' 
169 
170      END SELECT 
171
172     !    3.d.1 mbathy ( vertical k-levels of bathymetry )     
173
174     CALL crs_dom_bat
175     
176     !
177     CALL wrk_alloc(jpi, jpj, jpk, zfse3t, zfse3u, zfse3v, zfse3w )
178     !
179     zfse3t(:,:,:) = fse3t(:,:,:)
180     zfse3u(:,:,:) = fse3u(:,:,:)
181     zfse3v(:,:,:) = fse3v(:,:,:)
182     zfse3w(:,:,:) = fse3w(:,:,:)
183
184     !    3.d.2   Surfaces
185     CALL crs_dom_sfc( tmask, 'W', e1e2w_crs, e1e2w_msk, p_e1=e1t, p_e2=e2t    )
186     CALL crs_dom_sfc( umask, 'U', e2e3u_crs, e2e3u_msk, p_e2=e2u, p_e3=zfse3u )
187     CALL crs_dom_sfc( vmask, 'V', e1e3v_crs, e1e3v_msk, p_e1=e1v, p_e3=zfse3v )
188   
189     facsurfu(:,:,:) = umask_crs(:,:,:) * e2e3u_msk(:,:,:) / e2e3u_crs(:,:,:)
190     facsurfv(:,:,:) = vmask_crs(:,:,:) * e1e3v_msk(:,:,:) / e1e3v_crs(:,:,:)
191
192     !    3.d.3   Vertical scale factors
193     !
194   
195 
196     CALL crs_dom_e3( e1t, e2t, zfse3t, e1e2w_crs, 'T', tmask, e3t_crs, e3t_max_crs)
197     CALL crs_dom_e3( e1u, e2u, zfse3u, e2e3u_crs, 'U', umask, e3u_crs, e3u_max_crs)
198     CALL crs_dom_e3( e1v, e2v, zfse3v, e1e3v_crs, 'V', vmask, e3v_crs, e3v_max_crs)
199     CALL crs_dom_e3( e1t, e2t, zfse3w, e1e2w_crs, 'W', tmask, e3w_crs, e3w_max_crs)
200
201     ! Reset 0 to e3t_0 or e3w_0
202     DO jk = 1, jpk
203        DO ji = 1, jpi_crs
204           DO jj = 1, jpj_crs
205              IF( e3t_crs(ji,jj,jk) == 0._wp ) e3t_crs(ji,jj,jk) = e3t_0(jk)
206              IF( e3w_crs(ji,jj,jk) == 0._wp ) e3w_crs(ji,jj,jk) = e3w_0(jk)
207              IF( e3u_crs(ji,jj,jk) == 0._wp ) e3u_crs(ji,jj,jk) = e3t_0(jk)
208              IF( e3v_crs(ji,jj,jk) == 0._wp ) e3v_crs(ji,jj,jk) = e3t_0(jk)
209           ENDDO
210        ENDDO
211     ENDDO
212
213     !    3.d.3   Vertical depth (meters)
214     CALL crs_dom_ope( gdept, 'MAX', 'T', tmask, gdept_crs, p_e3=zfse3t ) 
215     CALL crs_dom_ope( gdepw, 'MAX', 'W', tmask, gdepw_crs, p_e3=zfse3w )
216
217
218     !---------------------------------------------------------
219     ! 4.  Coarse grid ocean volume and averaging weights
220     !---------------------------------------------------------
221     ! 4.a. Ocean volume or area unmasked and masked
222     CALL crs_dom_facvol( tmask, 'T', e1t, e2t, zfse3t, ocean_volume_crs_t, facvol_t )
223     !
224     bt_crs(:,:,:) = ocean_volume_crs_t(:,:,:) * facvol_t(:,:,:)
225     !
226     r1_bt_crs(:,:,:) = 0._wp 
227     WHERE( bt_crs /= 0._wp ) r1_bt_crs(:,:,:) = 1._wp / bt_crs(:,:,:)
228
229     CALL crs_dom_facvol( tmask, 'W', e1t, e2t, zfse3w, ocean_volume_crs_w, facvol_w )
230     !
231     !---------------------------------------------------------
232     ! 5.  Write out coarse meshmask  (see OPA_SRC/DOM/domwri.F90 for ideas later)
233     !---------------------------------------------------------
234
235     IF( nn_msh_crs > 0 ) THEN
236        CALL dom_grid_crs   ! Save the parent grid information  & Switch to coarse grid domain
237        CALL crs_dom_wri     
238        CALL dom_grid_glo   ! Return to parent grid domain
239     ENDIF
240     
241     !---------------------------------------------------------
242     ! 7. Finish and clean-up
243     !---------------------------------------------------------
244     CALL wrk_dealloc(jpi, jpj, jpk, zfse3t, zfse3u, zfse3v, zfse3w )
245
246
247   END SUBROUTINE crs_init
248   
249   !!======================================================================
250
251END MODULE crsini
Note: See TracBrowser for help on using the repository browser.