source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/CRS/crsini.F90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 5 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

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