source: NEMO/trunk/src/OCE/CRS/crsini.F90 @ 13237

Last change on this file since 13237 was 13237, checked in by smasson, 3 months ago

trunk: Mid-year merge, merge back KERNEL-06_techene_e3

  • Property svn:keywords set to Id
File size: 11.3 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   !!----------------------------------------------------------------------
10   !!  crs_init    :
11   !!----------------------------------------------------------------------
12   USE par_kind, ONLY: wp
13   USE par_oce                  ! For parameter jpi,jpj
14   USE dom_oce                  ! For parameters in par_oce
15   USE crs                      ! Coarse grid domain
16   USE phycst, ONLY: omega, rad ! physical constants
17   USE crsdom
18   USE crsdomwri
19   USE crslbclnk
20   !
21   USE iom
22   USE in_out_manager
23   USE lib_mpp
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   crs_init   ! called by nemogcm.F90 module
29
30   !! * Substitutions
31#  include "domzgr_substitute.h90"
32   !!----------------------------------------------------------------------
33   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
34   !! $Id$
35   !! Software governed by the CeCILL license (see ./LICENSE)
36   !!----------------------------------------------------------------------
37CONTAINS
38   
39   SUBROUTINE crs_init( Kmm )
40      !!-------------------------------------------------------------------
41      !!                     *** SUBROUTINE crs_init
42      !!  ** Purpose : Initialization of the grid coarsening module 
43      !!               1. Read namelist
44      !!               X2. MOVE TO crs_dom.F90 Set the domain definitions for coarse grid
45      !!                 a. Define the coarse grid starting/ending indices on parent grid
46      !!                    Here is where the T-pivot or F-pivot grids are discerned
47      !!                 b. TODO.  Include option for north-centric or equator-centric binning.
48      !!                 (centered only for odd reduction factors; even reduction bins bias equator to the south)
49      !!               3. Mask and mesh creation. => calls to crsfun
50      !!                  a. Use crsfun_mask to generate tmask,umask, vmask, fmask.
51      !!                  b. Use crsfun_coordinates to get coordinates
52      !!                  c. Use crsfun_UV to get horizontal scale factors
53      !!                  d. Use crsfun_TW to get initial vertical scale factors   
54      !!               4. Volumes and weights jes.... TODO. Updates for vvl? Where to do this? crsstp.F90?
55      !!                  a. Calculate initial coarse grid box volumes: pvol_T, pvol_W
56      !!                  b. Calculate initial coarse grid surface-averaging weights
57      !!                  c. Calculate initial coarse grid volume-averaging weights
58      !!                 
59      !!               X5. MOVE TO crs_dom_wri.F90 Using iom_rstput output the initial meshmask.
60      !!               ?. Another set of "masks" to generate
61      !!                  are the u- and v- surface areas for U- and V- area-weighted means.
62      !!                  Need to put this somewhere in section 3?
63      !! jes. What do to about the vvl?  GM.  could separate the weighting (denominator), so
64      !! output C*dA or C*dV as summation not mran, then do mean (division) at moment of output.
65      !! As is, crsfun takes into account vvl.   
66      !!      Talked about pre-setting the surface array to avoid IF/ENDIF and division.
67      !!      But have then to make that preset array here and elsewhere.
68      !!      that is called every timestep...
69      !!
70      !!               - Read in pertinent data ?
71      !!-------------------------------------------------------------------
72      INTEGER, INTENT(in) :: Kmm   ! time level index
73      !
74      INTEGER  :: ji,jj,jk      ! dummy indices
75      INTEGER  :: ierr                                ! allocation error status
76      INTEGER  ::   ios                 ! Local integer output status for namelist read
77      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, ze3u, ze3v, ze3w
78
79      NAMELIST/namcrs/ nn_factx, nn_facty, nn_binref, ln_msh_crs, nn_crs_kz, ln_crs_wn
80      !!----------------------------------------------------------------------
81      !
82     !---------------------------------------------------------
83     ! 1. Read Namelist file
84     !---------------------------------------------------------
85     !
86      READ  ( numnam_ref, namcrs, IOSTAT = ios, ERR = 901)
87901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namcrs in reference namelist' )
88      READ  ( numnam_cfg, namcrs, IOSTAT = ios, ERR = 902 )
89902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namcrs in configuration namelist' )
90      IF(lwm) WRITE ( numond, namcrs )
91
92     IF(lwp) THEN
93        WRITE(numout,*)
94        WRITE(numout,*) 'crs_init : Initializing the grid coarsening module'
95        WRITE(numout,*) '~~~~~~~~'
96        WRITE(numout,*) '   Namelist namcrs '
97        WRITE(numout,*) '      coarsening factor in i-direction      nn_factx   = ', nn_factx
98        WRITE(numout,*) '      coarsening factor in j-direction      nn_facty   = ', nn_facty
99        WRITE(numout,*) '      bin centering preference              nn_binref  = ', nn_binref
100        WRITE(numout,*) '      create a mesh file (=T)               ln_msh_crs = ', ln_msh_crs
101        WRITE(numout,*) '      type of Kz coarsening (0,1,2)         nn_crs_kz  = ', nn_crs_kz
102        WRITE(numout,*) '      ww coarsened or computed using hdiv  ln_crs_wn  = ', ln_crs_wn
103     ENDIF
104             
105     rfactx_r = 1. / nn_factx
106     rfacty_r = 1. / nn_facty
107
108     !---------------------------------------------------------
109     ! 2. Define Global Dimensions of the coarsened grid
110     !---------------------------------------------------------
111     CALL crs_dom_def     
112
113     !---------------------------------------------------------
114     ! 3. Mask and Mesh
115     !---------------------------------------------------------
116
117     !     Set up the masks and meshes     
118
119     !  3.a. Get the masks   
120 
121     CALL crs_dom_msk
122
123
124     !  3.b. Get the coordinates
125     !      Odd-numbered reduction factor, center coordinate on T-cell
126     !      Even-numbered reduction factor, center coordinate on U-,V- faces or f-corner.
127     !     
128     IF ( nresty /= 0 .AND. nrestx /= 0 ) THEN
129        CALL crs_dom_coordinates( gphit, glamt, 'T', gphit_crs, glamt_crs ) 
130        CALL crs_dom_coordinates( gphiu, glamu, 'U', gphiu_crs, glamu_crs )       
131        CALL crs_dom_coordinates( gphiv, glamv, '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( gphiu, glamu, 'T', gphit_crs, glamt_crs )
135        CALL crs_dom_coordinates( gphiu, glamu, 'U', gphiu_crs, glamu_crs )
136        CALL crs_dom_coordinates( gphif, glamf, '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( gphiv, glamv, 'T', gphit_crs, glamt_crs )
140        CALL crs_dom_coordinates( gphif, glamf, 'U', gphiu_crs, glamu_crs )
141        CALL crs_dom_coordinates( gphiv, glamv, 'V', gphiv_crs, glamv_crs )
142        CALL crs_dom_coordinates( gphif, glamf, 'F', gphif_crs, glamf_crs )
143     ELSE
144        CALL crs_dom_coordinates( gphif, glamf, 'T', gphit_crs, glamt_crs )
145        CALL crs_dom_coordinates( gphif, glamf, 'U', gphiu_crs, glamu_crs )
146        CALL crs_dom_coordinates( gphif, glamf, 'V', gphiv_crs, glamv_crs )
147        CALL crs_dom_coordinates( gphif, glamf, 'F', gphif_crs, glamf_crs )
148     ENDIF
149
150
151     !  3.c. Get the horizontal mesh
152
153     !      3.c.1 Horizontal scale factors
154
155     CALL crs_dom_hgr( e1t, e2t, 'T', e1t_crs, e2t_crs )
156     CALL crs_dom_hgr( e1u, e2u, 'U', e1u_crs, e2u_crs )
157     CALL crs_dom_hgr( e1v, e2v, 'V', e1v_crs, e2v_crs )
158     CALL crs_dom_hgr( e1f, e2f, 'F', e1f_crs, e2f_crs )
159
160     e1e2t_crs(:,:) = e1t_crs(:,:) * e2t_crs(:,:)
161     
162     
163     !      3.c.2 Coriolis factor 
164
165!!gm  Not sure CRS needs Coriolis parameter....
166!!gm  If needed, then update this to have Coriolis at both f- and t-points
167
168      ff_crs(:,:) = 2. * omega * SIN( rad * gphif_crs(:,:) )
169
170      CALL ctl_warn( 'crsini: CAUTION, CRS only designed for Coriolis defined on the sphere' ) 
171 
172
173     !    3.d.1 mbathy ( vertical k-levels of bathymetry )     
174
175     CALL crs_dom_bat
176     
177     !
178     DO jk = 1, jpk
179        ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
180        ze3u(:,:,jk) = e3u(:,:,jk,Kmm)
181        ze3v(:,:,jk) = e3v(:,:,jk,Kmm)
182        ze3w(:,:,jk) = e3w(:,:,jk,Kmm)
183     END DO 
184
185     !    3.d.2   Surfaces
186     CALL crs_dom_sfc( tmask, 'W', e1e2w_crs, e1e2w_msk, p_e1=e1t, p_e2=e2t  )
187     CALL crs_dom_sfc( umask, 'U', e2e3u_crs, e2e3u_msk, p_e2=e2u, p_e3=ze3u )
188     CALL crs_dom_sfc( vmask, 'V', e1e3v_crs, e1e3v_msk, p_e1=e1v, p_e3=ze3v )
189   
190     facsurfu(:,:,:) = umask_crs(:,:,:) * e2e3u_msk(:,:,:) / e2e3u_crs(:,:,:)
191     facsurfv(:,:,:) = vmask_crs(:,:,:) * e1e3v_msk(:,:,:) / e1e3v_crs(:,:,:)
192
193     !    3.d.3   Vertical scale factors
194     !
195     CALL crs_dom_e3( e1t, e2t, ze3t, e1e2w_crs, 'T', tmask, e3t_crs, e3t_max_crs)
196     CALL crs_dom_e3( e1u, e2u, ze3u, e2e3u_crs, 'U', umask, e3u_crs, e3u_max_crs)
197     CALL crs_dom_e3( e1v, e2v, ze3v, e1e3v_crs, 'V', vmask, e3v_crs, e3v_max_crs)
198     CALL crs_dom_e3( e1t, e2t, ze3w, e1e2w_crs, 'W', tmask, e3w_crs, e3w_max_crs)
199
200     ! Replace 0 by e3t_0 or e3w_0
201     DO jk = 1, jpk
202        DO ji = 1, jpi_crs
203           DO jj = 1, jpj_crs
204              IF( e3t_crs(ji,jj,jk) == 0._wp ) e3t_crs(ji,jj,jk) = e3t_1d(jk)
205              IF( e3w_crs(ji,jj,jk) == 0._wp ) e3w_crs(ji,jj,jk) = e3w_1d(jk)
206              IF( e3u_crs(ji,jj,jk) == 0._wp ) e3u_crs(ji,jj,jk) = e3t_1d(jk)
207              IF( e3v_crs(ji,jj,jk) == 0._wp ) e3v_crs(ji,jj,jk) = e3t_1d(jk)
208           ENDDO
209        ENDDO
210     ENDDO
211
212     !    3.d.3   Vertical depth (meters)
213     CALL crs_dom_ope( gdept_0, 'MAX', 'T', tmask, gdept_crs, p_e3=ze3t, psgn=1.0_wp ) 
214     CALL crs_dom_ope( gdepw_0, 'MAX', 'W', tmask, gdepw_crs, p_e3=ze3w, psgn=1.0_wp )
215
216
217     !---------------------------------------------------------
218     ! 4.  Coarse grid ocean volume and averaging weights
219     !---------------------------------------------------------
220     ! 4.a. Ocean volume or area unmasked and masked
221     CALL crs_dom_facvol( tmask, 'T', e1t, e2t, ze3t, ocean_volume_crs_t, facvol_t )
222     !
223     bt_crs(:,:,:) = ocean_volume_crs_t(:,:,:) * facvol_t(:,:,:)
224     !
225     r1_bt_crs(:,:,:) = 0._wp 
226     WHERE( bt_crs /= 0._wp ) r1_bt_crs(:,:,:) = 1._wp / bt_crs(:,:,:)
227
228     CALL crs_dom_facvol( tmask, 'W', e1t, e2t, ze3w, ocean_volume_crs_w, facvol_w )
229     !
230     !---------------------------------------------------------
231     ! 5.  Write out coarse meshmask  (see OCE/DOM/domwri.F90 for ideas later)
232     !---------------------------------------------------------
233
234     IF( ln_msh_crs ) THEN
235        CALL dom_grid_crs   ! Save the parent grid information  & Switch to coarse grid domain
236        CALL crs_dom_wri     
237        CALL dom_grid_glo   ! Return to parent grid domain
238     ENDIF
239     
240      !---------------------------------------------------------
241      ! 7. Finish and clean-up
242      !---------------------------------------------------------
243      !
244   END SUBROUTINE crs_init
245   
246   !!======================================================================
247END MODULE crsini
Note: See TracBrowser for help on using the repository browser.