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/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/CRS – NEMO

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/CRS/crsini.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago

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