1 | MODULE crsini |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE crsini *** |
---|
4 | !! Manage the grid coarsening module initialization |
---|
5 | !!====================================================================== |
---|
6 | !! History 2012-05 (J. Simeon, G. Madec, C. Ethe) 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_dom ! 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 crs |
---|
18 | USE crsdomwri |
---|
19 | USE crslbclnk |
---|
20 | |
---|
21 | IMPLICIT NONE |
---|
22 | PRIVATE |
---|
23 | |
---|
24 | PUBLIC crs_init |
---|
25 | |
---|
26 | !! * Substitutions |
---|
27 | # include "domzgr_substitute.h90" |
---|
28 | |
---|
29 | CONTAINS |
---|
30 | |
---|
31 | SUBROUTINE crs_init |
---|
32 | !!------------------------------------------------------------------- |
---|
33 | !! *** SUBROUTINE crs_init |
---|
34 | !! ** Purpose : Initialization of the grid coarsening module |
---|
35 | !! 1. Read namelist |
---|
36 | !! X2. MOVE TO crs_dom.F90 Set the domain definitions for coarse grid |
---|
37 | !! a. Define the coarse grid starting/ending indices on parent grid |
---|
38 | !! Here is where the T-pivot or F-pivot grids are discerned |
---|
39 | !! b. TODO. Include option for north-centric or equator-centric binning. |
---|
40 | !! (centered only for odd reduction factors; even reduction bins bias equator to the south) |
---|
41 | !! 3. Mask and mesh creation. => calls to crsfun |
---|
42 | !! a. Use crsfun_mask to generate tmask,umask, vmask, fmask. |
---|
43 | !! b. Use crsfun_coordinates to get coordinates |
---|
44 | !! c. Use crsfun_UV to get horizontal scale factors |
---|
45 | !! d. Use crsfun_TW to get initial vertical scale factors |
---|
46 | !! 4. Volumes and weights jes.... TODO. Updates for vvl? Where to do this? crsstp.F90? |
---|
47 | !! a. Calculate initial coarse grid box volumes: pvol_T, pvol_W |
---|
48 | !! b. Calculate initial coarse grid surface-averaging weights |
---|
49 | !! c. Calculate initial coarse grid volume-averaging weights |
---|
50 | !! |
---|
51 | !! X5. MOVE TO crs_dom_wri.F90 Using iom_rstput output the initial meshmask. |
---|
52 | !! ?. Another set of "masks" to generate |
---|
53 | !! are the u- and v- surface areas for U- and V- area-weighted means. |
---|
54 | !! Need to put this somewhere in section 3? |
---|
55 | !! jes. What do to about the vvl? GM. could separate the weighting (denominator), so |
---|
56 | !! output C*dA or C*dV as summation not mran, then do mean (division) at moment of output. |
---|
57 | !! As is, crsfun takes into account vvl. |
---|
58 | !! Talked about pre-setting the surface array to avoid IF/ENDIFS and division. |
---|
59 | !! But have then to make that preset array here and elsewhere. |
---|
60 | !! that is called every timestep... |
---|
61 | !! |
---|
62 | !! - Read in pertinent data ? |
---|
63 | !!------------------------------------------------------------------- |
---|
64 | !! Local variables |
---|
65 | INTEGER :: ji,jj,jk,ijjgloT,ijis,ijie,ijjs,ijje ! dummy indices |
---|
66 | INTEGER :: ierr ! allocation error status |
---|
67 | REAL(wp) :: zrestx, zresty ! for determining odd or even reduction factor |
---|
68 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zmbk |
---|
69 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zfse3t, zfse3u, zfse3v, zfse3f |
---|
70 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zfse3w, zfse3t_n, zfse3t_b |
---|
71 | LOGICAL :: llok |
---|
72 | |
---|
73 | NAMELIST/namcrs/ nn_factx, nn_facty, cn_binref, nn_fcrs, nn_msh_crs, cn_ocerstcrs |
---|
74 | |
---|
75 | !!---------------------------------------------------------------------- |
---|
76 | ! |
---|
77 | IF( nn_timing == 1 ) CALL timing_start('crs_init') |
---|
78 | ! |
---|
79 | IF(lwp) THEN |
---|
80 | WRITE(numout,*) |
---|
81 | WRITE(numout,*) 'crs_init : Initializing the grid coarsening module ' |
---|
82 | ENDIF |
---|
83 | |
---|
84 | !--------------------------------------------------------- |
---|
85 | ! 1. Read Namelist file |
---|
86 | !--------------------------------------------------------- |
---|
87 | ! |
---|
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,*) ' nn_factx = ', nn_factx |
---|
93 | WRITE(numout,*) ' nn_facty = ', nn_facty |
---|
94 | WRITE(numout,*) ' cn_binref = ', cn_binref |
---|
95 | WRITE(numout,*) ' nn_fcrs = ', nn_fcrs |
---|
96 | WRITE(numout,*) ' nn_msh_crs = ', nn_msh_crs |
---|
97 | ENDIF |
---|
98 | |
---|
99 | rfactx_r = 1./nn_factx |
---|
100 | rfacty_r = 1./nn_facty |
---|
101 | |
---|
102 | !--------------------------------------------------------- |
---|
103 | ! 2. Define Global Dimensions of the coarsened grid |
---|
104 | !--------------------------------------------------------- |
---|
105 | ! 2.a. Define global domain indices |
---|
106 | jpiglo_crs = INT( (jpiglo - 2) / nn_factx ) + 2 |
---|
107 | jpjglo_crs = INT( (jpjglo - 2) / nn_facty ) + 2 ! the -2 removes j=1, j=jpj |
---|
108 | jpiglo_crsm1 = jpiglo_crs - 1 |
---|
109 | jpjglo_crsm1 = jpjglo_crs - 1 |
---|
110 | jpkm1 = jpk - 1 |
---|
111 | |
---|
112 | ! 2.b. Define local domain indices |
---|
113 | jpi_crs = ( jpiglo_crs-2*jpreci + (jpni-1) ) / jpni + 2*jpreci |
---|
114 | jpj_crs = ( jpjglo_crs-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj |
---|
115 | jpi_crsm1 = jpi_crs - 1 |
---|
116 | jpj_crsm1 = jpj_crs - 1 |
---|
117 | |
---|
118 | nperio_crs = jperio |
---|
119 | npolj_crs = npolj |
---|
120 | |
---|
121 | IF ( jpnij == 1 ) THEN |
---|
122 | jpnij_crs = jpnij |
---|
123 | narea_crs = narea |
---|
124 | nimpp_crs = nimpp |
---|
125 | njmpp_crs = njmpp |
---|
126 | ELSE |
---|
127 | WRITE(numout,*) 'crsini.F90. mpp not supported... Stopping' |
---|
128 | STOP |
---|
129 | ENDIF |
---|
130 | |
---|
131 | nlcj_crs = jpj_crs |
---|
132 | nlci_crs = jpi_crs |
---|
133 | nldi_crs = 1 |
---|
134 | nlei_crs = jpi_crs |
---|
135 | nlej_crs = jpj_crs |
---|
136 | nldj_crs = 1 |
---|
137 | |
---|
138 | IF (lwp) THEN |
---|
139 | WRITE(numout,*) |
---|
140 | WRITE(numout,*) 'crs_init : coarse grid dimensions' |
---|
141 | WRITE(numout,*) '~~~~~~~ coarse domain global j-dimension jpjglo_crs = ', jpjglo_crs |
---|
142 | WRITE(numout,*) '~~~~~~~ coarse domain global i-dimension jpiglo_crs = ', jpiglo_crs |
---|
143 | WRITE(numout,*) '~~~~~~~ coarse domain local i-dimension jpi_crs = ', jpi_crs |
---|
144 | WRITE(numout,*) '~~~~~~~ coarse domain local j-dimension jpj_crs = ', jpj_crs |
---|
145 | ENDIF |
---|
146 | |
---|
147 | |
---|
148 | mxbinctr = INT( nn_factx * 0.5 ) |
---|
149 | mybinctr = INT( nn_facty * 0.5 ) |
---|
150 | |
---|
151 | zrestx = MOD( nn_factx, 2 ) ! check if even- or odd- numbered reduction factor |
---|
152 | zresty = MOD( nn_facty, 2 ) |
---|
153 | |
---|
154 | IF ( zrestx == 0 ) THEN |
---|
155 | mxbinctr = mxbinctr - 1 |
---|
156 | ENDIF |
---|
157 | |
---|
158 | IF ( zresty == 0 ) THEN |
---|
159 | mybinctr = mybinctr - 1 |
---|
160 | IF ( jperio == 3 .OR. jperio == 4 ) nperio_crs = jperio + 2 |
---|
161 | IF ( jperio == 5 .OR. jperio == 6 ) nperio_crs = jperio - 2 |
---|
162 | |
---|
163 | IF ( npolj == 3 ) npolj_crs = 5 |
---|
164 | IF ( npolj == 5 ) npolj_crs = 3 |
---|
165 | ENDIF |
---|
166 | |
---|
167 | rfactxy = nn_factx * nn_facty |
---|
168 | |
---|
169 | |
---|
170 | !jes. TODO Need to deallocate these if ln_crs = F |
---|
171 | ierr = crs_dom_alloc() ! allocate most coarse grid arrays |
---|
172 | |
---|
173 | ! jes. TODO. Add the next two lines when mpp is done |
---|
174 | ! IF( lk_mpp ) CALL mpp_sum( ierr ) |
---|
175 | ! IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'nemo_alloc : unable to allocate standard ocean arrays' ) |
---|
176 | |
---|
177 | |
---|
178 | ! 2.c. Set up bins for coarse grid, horizontal only. |
---|
179 | |
---|
180 | mis_crs(:) = 0; mie_crs(:) = 0 |
---|
181 | mjs_crs(:) = 0; mje_crs(:) = 0 |
---|
182 | |
---|
183 | SELECT CASE ( cn_binref ) |
---|
184 | |
---|
185 | CASE ( 'NORTH' ) |
---|
186 | |
---|
187 | SELECT CASE ( nperio ) |
---|
188 | |
---|
189 | CASE ( 2 ) |
---|
190 | WRITE(numout,*) 'crs_init, jperio=2 not supported' |
---|
191 | |
---|
192 | CASE ( 3, 4 ) ! T-Pivot at North Fold |
---|
193 | |
---|
194 | DO ji = 2, jpiglo_crsm1 |
---|
195 | !cc ijie = (ji*nn_factx)-nn_factx+1 |
---|
196 | ijie = (ji*nn_factx)-nn_factx !cc |
---|
197 | ijis = ijie-nn_factx+1 |
---|
198 | |
---|
199 | IF ( ji == jpiglo_crsm1 ) THEN |
---|
200 | IF ( ((jpiglo-1)-ijie) <= nn_factx ) ijie = jpiglo-2 ! ijie = jpiglo-1 !cc |
---|
201 | ENDIF |
---|
202 | |
---|
203 | ! Handle first the northernmost bin |
---|
204 | IF ( nn_facty == 2 ) THEN |
---|
205 | ijjgloT=jpjglo-1 |
---|
206 | ELSE |
---|
207 | ijjgloT=jpjglo |
---|
208 | ENDIF |
---|
209 | |
---|
210 | DO jj = 2, jpjglo_crsm1 |
---|
211 | ! cc ijje = ijjgloT-nn_facty*(jj-2) |
---|
212 | ijje = ijjgloT-nn_facty*(jj-2) - 1 |
---|
213 | ijjs = ijje-nn_facty+1 |
---|
214 | |
---|
215 | IF ( ijjs <= nn_facty ) ijjs = 2 |
---|
216 | |
---|
217 | mis_crs(ji) = ijis |
---|
218 | mie_crs(ji) = ijie |
---|
219 | mjs_crs(jpjglo_crs-jj+1) = ijjs |
---|
220 | mje_crs(jpjglo_crs-jj+1) = ijje |
---|
221 | |
---|
222 | ENDDO |
---|
223 | ENDDO |
---|
224 | |
---|
225 | CASE ( 5, 6 ) ! F-pivot at North Fold |
---|
226 | |
---|
227 | DO ji = 2, jpiglo_crsm1 |
---|
228 | ijie = (ji*nn_factx)-nn_factx+1 |
---|
229 | ijis = ijie-nn_factx+1 |
---|
230 | |
---|
231 | IF ( ji == jpiglo_crsm1 ) THEN |
---|
232 | IF ( ((jpiglo-1)-ijie) <= nn_factx ) ijie = jpiglo-1 |
---|
233 | ENDIF |
---|
234 | |
---|
235 | ! Treat the northernmost bin separately. |
---|
236 | jj = 2 |
---|
237 | ijje = jpjglo-nn_facty*(jj-2) |
---|
238 | IF ( nn_facty == 3 ) THEN |
---|
239 | ijjs=ijje-1 |
---|
240 | ELSE |
---|
241 | ijjs=ijje-nn_facty+1 |
---|
242 | ENDIF |
---|
243 | |
---|
244 | mis_crs(ji) = ijis |
---|
245 | mie_crs(ji) = ijie |
---|
246 | mjs_crs(jpjglo_crs-jj+1) = ijjs |
---|
247 | mje_crs(jpjglo_crs-jj+1) = ijje |
---|
248 | |
---|
249 | ! Now bin the rest, any remainder at the south is lumped in the southern bin |
---|
250 | DO jj = 3, jpjglo_crsm1 |
---|
251 | |
---|
252 | ijje = jpjglo-nn_facty*(jj-2) |
---|
253 | ijjs = ijje-nn_facty+1 |
---|
254 | |
---|
255 | IF ( ijjs <= nn_facty ) ijjs = 2 |
---|
256 | |
---|
257 | mis_crs(ji) = ijis |
---|
258 | mie_crs(ji) = ijie |
---|
259 | mjs_crs(jpjglo_crs-jj+1) = ijjs |
---|
260 | mje_crs(jpjglo_crs-jj+1) = ijje |
---|
261 | ENDDO |
---|
262 | ENDDO |
---|
263 | |
---|
264 | CASE DEFAULT |
---|
265 | WRITE(numout,*) 'crs_init. Only jperio = 3, 4, 5, 6 supported' |
---|
266 | |
---|
267 | END SELECT |
---|
268 | |
---|
269 | CASE ( 'EQUAT' ) |
---|
270 | WRITE(numout,*) 'crs_init. Equator-centered bins option not yet available' |
---|
271 | |
---|
272 | END SELECT |
---|
273 | |
---|
274 | ! Pad the boundaries, do not know if it is necessary |
---|
275 | mis_crs(1) = 1 ; mis_crs(jpiglo_crs) = mie_crs(jpiglo_crs - 1) + 1 !cc |
---|
276 | mie_crs(1) = nn_factx ; mie_crs(jpiglo_crs) = jpiglo !cc |
---|
277 | mjs_crs(1) = 1 ; mjs_crs(jpjglo_crs) = mje_crs(jpjglo_crs - 1) + 1 |
---|
278 | mje_crs(1) = mjs_crs(2)-1; mje_crs(jpjglo_crs) = jpjglo |
---|
279 | |
---|
280 | ! WRITE(numout,*) 'crs_init. coarse grid bounds on parent grid' |
---|
281 | ! WRITE(numout,'(1x,a,62(1x,i3),/)') 'mis_crs=', mis_crs |
---|
282 | ! WRITE(numout,'(1x,a,62(1x,i3),/)') 'mie_crs=', mie_crs |
---|
283 | ! WRITE(numout,'(1x,a,51(1x,i3),/)') 'mjs_crs=', mjs_crs |
---|
284 | ! WRITE(numout,'(1x,a,51(1x,i3),/)') 'mje_crs=', mje_crs |
---|
285 | |
---|
286 | |
---|
287 | !--------------------------------------------------------- |
---|
288 | ! 3. Mask and Mesh |
---|
289 | !--------------------------------------------------------- |
---|
290 | |
---|
291 | ! Set up the masks and meshes |
---|
292 | |
---|
293 | ! 3.a. Get the masks |
---|
294 | CALL crsfun( p_ptmask=tmask, cd_type='T', p_cmask=tmask_crs ) |
---|
295 | CALL crsfun( p_ptmask=umask, cd_type='U', p_cmask=umask_crs ) |
---|
296 | CALL crsfun( p_ptmask=vmask, cd_type='V', p_cmask=vmask_crs ) |
---|
297 | CALL crsfun( p_ptmask=fmask, cd_type='F', p_cmask=fmask_crs ) |
---|
298 | |
---|
299 | |
---|
300 | ! CALL crsfun( p_ptmask=tmask, cd_type='T', p_pmask2d=rnfmsk, p_cmask2d=rnfmsk_crs ) |
---|
301 | ! rnfmsk_crs(:,:) = rnfmsk_crs(:,:) * tmask_crs(:,:,1) |
---|
302 | |
---|
303 | WRITE(numout,*) 'crsini. Finished masks' |
---|
304 | |
---|
305 | ! 3.b. Get the coordinates |
---|
306 | ! Odd-numbered reduction factor, center coordinate on T-cell |
---|
307 | ! Even-numbered reduction factor, center coordinate on U-,V- faces or f-corner. |
---|
308 | ! |
---|
309 | IF ( zresty /= 0 .AND. zrestx /= 0 ) THEN |
---|
310 | |
---|
311 | CALL crsfun( gphit, glamt, 'T', gphit_crs, glamt_crs ) |
---|
312 | WRITE(numout,*) 'crsini. gphit_crs(15,15)', gphit_crs(15,15) |
---|
313 | WRITE(numout,*) 'crsini. glamt_crs(15,15)', glamt_crs(15,15) |
---|
314 | |
---|
315 | WRITE(numout,*) 'crsini. count 1' |
---|
316 | |
---|
317 | CALL crsfun( gphiu, glamu, 'U', gphiu_crs, glamu_crs ) !cc |
---|
318 | WRITE(numout,*) 'crsini. gphiu_crs(15,15)', gphiu_crs(15,15) !cc |
---|
319 | WRITE(numout,*) 'crsini. glamu_crs(15,15)', glamu_crs(15,15) !cc |
---|
320 | WRITE(numout,*) 'crsini. count 2' |
---|
321 | |
---|
322 | CALL crsfun( p_pgphi=gphiv, p_pglam=glamv, cd_type='V', p_cgphi=gphiv_crs, p_cglam=glamv_crs ) !cc |
---|
323 | WRITE(numout,*) 'crsini. gphiv_crs(15,15)', gphiv_crs(15,15) !cc |
---|
324 | WRITE(numout,*) 'crsini. glamv_crs(15,15)', glamv_crs(15,15) !cc |
---|
325 | |
---|
326 | WRITE(numout,*) 'crsini. count 3' |
---|
327 | CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='F', p_cgphi=gphif_crs, p_cglam=glamf_crs ) !cc |
---|
328 | WRITE(numout,*) 'crsini. gphif_crs(15,15)', gphif_crs(15,15) !cc |
---|
329 | WRITE(numout,*) 'crsini. glamf_crs(15,15)', glamf_crs(15,15) !cc |
---|
330 | |
---|
331 | WRITE(numout,*) 'crsini. count 4' |
---|
332 | ELSEIF ( zresty /= 0 .AND. zrestx == 0 ) THEN |
---|
333 | CALL crsfun( p_pgphi=gphiu, p_pglam=glamu, cd_type='T', p_cgphi=gphit_crs, p_cglam=glamt_crs ) |
---|
334 | CALL crsfun( p_pgphi=gphiu, p_pglam=glamu, cd_type='U', p_cgphi=gphiu_crs, p_cglam=glamu_crs ) |
---|
335 | CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='V', p_cgphi=gphiv_crs, p_cglam=glamv_crs ) |
---|
336 | CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='F', p_cgphi=gphif_crs, p_cglam=glamf_crs ) |
---|
337 | ELSEIF ( zresty == 0 .AND. zrestx /= 0 ) THEN |
---|
338 | CALL crsfun( p_pgphi=gphiv, p_pglam=glamv, cd_type='T', p_cgphi=gphit_crs, p_cglam=glamt_crs ) |
---|
339 | CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='U', p_cgphi=gphiu_crs, p_cglam=glamu_crs ) |
---|
340 | CALL crsfun( p_pgphi=gphiv, p_pglam=glamv, cd_type='V', p_cgphi=gphiv_crs, p_cglam=glamv_crs ) |
---|
341 | CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='F', p_cgphi=gphif_crs, p_cglam=glamf_crs ) |
---|
342 | ELSE |
---|
343 | CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='T', p_cgphi=gphit_crs, p_cglam=glamt_crs ) |
---|
344 | CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='U', p_cgphi=gphiu_crs, p_cglam=glamu_crs ) |
---|
345 | CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='V', p_cgphi=gphiv_crs, p_cglam=glamv_crs ) |
---|
346 | CALL crsfun( p_pgphi=gphif, p_pglam=glamf, cd_type='F', p_cgphi=gphif_crs, p_cglam=glamf_crs ) |
---|
347 | ENDIF |
---|
348 | |
---|
349 | WRITE(numout,*) 'crsini. Finished coordinates' |
---|
350 | |
---|
351 | ! 3.c. Get the horizontal mesh |
---|
352 | |
---|
353 | ! 3.c.1 Horizontal scale factors |
---|
354 | ! CALL crsfun( cd_type='T', cd_op='SUM', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_cfield2d_1=e1t_crs, p_cfield2d_2=e2t_crs ) |
---|
355 | ! CALL crsfun( cd_type='U', cd_op='SUM', p_pmask=umask, p_e1=e1u, p_e2=e2u, p_cfield2d_1=e1u_crs, p_cfield2d_2=e2u_crs ) |
---|
356 | ! CALL crsfun( cd_type='V', cd_op='SUM', p_pmask=vmask, p_e1=e1v, p_e2=e2v, p_cfield2d_1=e1v_crs, p_cfield2d_2=e2v_crs ) |
---|
357 | ! CALL crsfun( cd_type='F', cd_op='SUM', p_pmask=fmask, p_e1=e1f, p_e2=e2f, p_cfield2d_1=e1f_crs, p_cfield2d_2=e2f_crs ) |
---|
358 | CALL crsfun( cd_type='T', cd_op='POS', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_cfield2d_1=e1t_crs, p_cfield2d_2=e2t_crs ) |
---|
359 | CALL crsfun( cd_type='U', cd_op='POS', p_pmask=umask, p_e1=e1u, p_e2=e2u, p_cfield2d_1=e1u_crs, p_cfield2d_2=e2u_crs ) |
---|
360 | CALL crsfun( cd_type='V', cd_op='POS', p_pmask=vmask, p_e1=e1v, p_e2=e2v, p_cfield2d_1=e1v_crs, p_cfield2d_2=e2v_crs ) |
---|
361 | CALL crsfun( cd_type='F', cd_op='POS', p_pmask=fmask, p_e1=e1f, p_e2=e2f, p_cfield2d_1=e1f_crs, p_cfield2d_2=e2f_crs ) |
---|
362 | |
---|
363 | e1e2t_crs(:,:) = e1t_crs(:,:) * e2t_crs(:,:) |
---|
364 | |
---|
365 | WRITE(numout,*) 'crsini. Finished horizontal scale factors' |
---|
366 | |
---|
367 | ! 3.c.2 Coriolis factor |
---|
368 | |
---|
369 | SELECT CASE( jphgr_msh ) ! type of horizontal mesh |
---|
370 | |
---|
371 | CASE ( 0, 1, 4 ) ! mesh on the sphere |
---|
372 | |
---|
373 | ff_crs(:,:) = 2. * omega * SIN( rad * gphif_crs(:,:) ) |
---|
374 | |
---|
375 | CASE DEFAULT |
---|
376 | |
---|
377 | WRITE(numout,*) 'crsini.F90. crs_init. Only jphgr_msh = 0, 1 or 4 supported' |
---|
378 | |
---|
379 | END SELECT |
---|
380 | |
---|
381 | |
---|
382 | ! 3.d. Get the initial vertical mesh |
---|
383 | ! nav_lon(jpi,jpj), nav_lat(jpi,jpj) |
---|
384 | ! nav_lev(jpk), e3t_0(jpk), e3w_0(jpk), gdep[t,w]_0(jpk) -> these stay constant |
---|
385 | ! 2D: mbathy (k-levels) |
---|
386 | ! 3D: fse3[t,u,v,w,f] (scale factors), gdep[t,u,v,w] (depth in meters) |
---|
387 | |
---|
388 | ! jes. TODO. Could probably make optional e1e2t in crsfun_TW... |
---|
389 | |
---|
390 | ! 3.d.1 mbathy ( vertical k-levels of bathymetry ) |
---|
391 | |
---|
392 | mbathy_crs(:,:) = jpkm1 |
---|
393 | mbkt_crs(:,:) = 1 |
---|
394 | mbku_crs(:,:) = 1 |
---|
395 | mbkv_crs(:,:) = 1 |
---|
396 | |
---|
397 | |
---|
398 | DO jj = 1, jpj_crs |
---|
399 | DO ji = 1, jpi_crs |
---|
400 | jk = 0 |
---|
401 | DO WHILE( tmask_crs(ji,jj,jk+1) > 0.) |
---|
402 | jk = jk + 1 |
---|
403 | ENDDO |
---|
404 | mbathy_crs(ji,jj) = float( jk ) |
---|
405 | ENDDO |
---|
406 | ENDDO |
---|
407 | |
---|
408 | ALLOCATE( zmbk(jpi_crs,jpj_crs) ) |
---|
409 | |
---|
410 | zmbk(:,:) = 0.0 |
---|
411 | zmbk(:,:) = REAL( mbathy_crs(:,:), wp ) ; CALL crs_lbc_lnk('T',1.0,zmbk) ; mbathy_crs(:,:) = INT( zmbk(:,:) ) |
---|
412 | |
---|
413 | |
---|
414 | ! |
---|
415 | IF(lwp) WRITE(numout,*) |
---|
416 | IF(lwp) WRITE(numout,*) ' crsini : mbkt is ocean bottom k-index of T-, U-, V- and W-levels ' |
---|
417 | IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' |
---|
418 | ! |
---|
419 | mbkt_crs(:,:) = MAX( mbathy_crs(:,:) , 1 ) ! bottom k-index of T-level (=1 over land) |
---|
420 | ! ! bottom k-index of W-level = mbkt+1 |
---|
421 | |
---|
422 | DO jj = 1, jpj_crsm1 ! bottom k-index of u- (v-) level |
---|
423 | DO ji = 1, jpi_crsm1 |
---|
424 | mbku_crs(ji,jj) = MIN( mbkt_crs(ji+1,jj ) , mbkt_crs(ji,jj) ) |
---|
425 | mbkv_crs(ji,jj) = MIN( mbkt_crs(ji ,jj+1) , mbkt_crs(ji,jj) ) |
---|
426 | END DO |
---|
427 | END DO |
---|
428 | |
---|
429 | WRITE(numout,*) 'crsini. Set mbku, mkbv' |
---|
430 | |
---|
431 | ! convert into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk |
---|
432 | zmbk(:,:) = 1.e0 |
---|
433 | zmbk(:,:) = REAL( mbku_crs(:,:), wp ) ; CALL crs_lbc_lnk('U',1.0,zmbk) ; mbku_crs (:,:) = MAX( INT( zmbk(:,:) ), 1 ) |
---|
434 | zmbk(:,:) = REAL( mbkv_crs(:,:), wp ) ; CALL crs_lbc_lnk('V',1.0,zmbk) ; mbkv_crs (:,:) = MAX( INT( zmbk(:,:) ), 1 ) |
---|
435 | ! |
---|
436 | |
---|
437 | WRITE(numout,*) 'crs_init = finished section 3d.1 jpi=', jpi, 'jpj=',jpj, ' jpk=', jpk |
---|
438 | |
---|
439 | ! 3.d.2 Vertical scale factors |
---|
440 | |
---|
441 | ALLOCATE( zfse3t(jpi,jpj,jpk), zfse3u(jpi,jpj,jpk), zfse3v(jpi,jpj,jpk), zfse3f(jpi,jpj,jpk), & |
---|
442 | & zfse3w(jpi,jpj,jpk), zfse3t_n(jpi,jpj,jpk), zfse3t_b(jpi,jpj,jpk) ) |
---|
443 | zfse3t(:,:,:) = fse3t(:,:,:) |
---|
444 | zfse3u(:,:,:) = fse3u(:,:,:) |
---|
445 | zfse3v(:,:,:) = fse3v(:,:,:) |
---|
446 | zfse3f(:,:,:) = fse3f(:,:,:) |
---|
447 | zfse3w(:,:,:) = fse3w(:,:,:) |
---|
448 | |
---|
449 | |
---|
450 | |
---|
451 | !CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='MAX', p_cmask=tmask_crs, p_ptmask=tmask, p_pfield3d_1=zfse3t, p_cfield3d=e3t_crs ) |
---|
452 | !CALL crsfun( p_e1e2t=e1e2t, cd_type='W', cd_op='MAX', p_cmask=tmask_crs, p_ptmask=tmask, p_pfield3d_1=zfse3w, p_cfield3d=e3w_crs ) |
---|
453 | !CALL crsfun( p_e1e2t=e1e2t, cd_type='U', cd_op='MIN', p_cmask=umask_crs, p_ptmask=umask, p_pfield3d_1=zfse3u, p_cfield3d=e3u_crs ) |
---|
454 | !CALL crsfun( p_e1e2t=e1e2t, cd_type='V', cd_op='MIN', p_cmask=vmask_crs, p_ptmask=vmask, p_pfield3d_1=zfse3v, p_cfield3d=e3v_crs ) |
---|
455 | !CALL crsfun( p_e1e2t=e1e2t, cd_type='F', cd_op='MIN', p_cmask=fmask_crs, p_ptmask=fmask, p_pfield3d_1=zfse3f, p_cfield3d=e3f_crs ) |
---|
456 | CALL crs_e3_max( p_e3=zfse3t, cd_type='T', p_mask=tmask, p_e3_crs=e3t_crs) |
---|
457 | CALL crs_e3_max( p_e3=zfse3w, cd_type='W', p_mask=tmask, p_e3_crs=e3w_crs) |
---|
458 | |
---|
459 | ! Reset 0 to e3t_0 or e3w_0 |
---|
460 | DO jk = 1, jpk |
---|
461 | DO ji = 1, jpi_crs |
---|
462 | DO jj = 1, jpj_crs |
---|
463 | IF ( e3t_crs(ji,jj,jk) == 0 ) e3t_crs(ji,jj,jk) = e3t_0(jk) |
---|
464 | IF ( e3w_crs(ji,jj,jk) == 0 ) e3w_crs(ji,jj,jk) = e3w_0(jk) |
---|
465 | IF ( e3u_crs(ji,jj,jk) == 0 ) e3u_crs(ji,jj,jk) = e3t_0(jk) |
---|
466 | IF ( e3v_crs(ji,jj,jk) == 0 ) e3v_crs(ji,jj,jk) = e3t_0(jk) |
---|
467 | IF ( e3f_crs(ji,jj,jk) == 0 ) e3f_crs(ji,jj,jk) = e3t_0(jk) |
---|
468 | ENDDO |
---|
469 | ENDDO |
---|
470 | ENDDO |
---|
471 | |
---|
472 | ! 3.d.3 Vertical depth (meters) |
---|
473 | |
---|
474 | CALL crsfun( p_e1e2t=e1e2t, cd_type='T', cd_op='MAX', p_cmask=tmask_crs, p_ptmask=tmask, & |
---|
475 | & p_pfield3d_1=gdept, p_cfield3d=gdept_crs ) |
---|
476 | CALL crsfun( p_e1e2t=e1e2t, cd_type='W', cd_op='MAX', p_cmask=tmask_crs, p_ptmask=tmask, & |
---|
477 | & p_pfield3d_1=gdepw, p_cfield3d=gdepw_crs ) |
---|
478 | |
---|
479 | ! 3.d.4 Surfaces |
---|
480 | |
---|
481 | CALL crs_surf(p_e1=e1t, p_e2=e2t ,p_e3=zfse3w, cd_type='W', p_mask=tmask, surf_crs=e1e2w, surf_msk_crs=e1e2w_msk) |
---|
482 | CALL crs_surf(p_e1=e1u, p_e2=e2u ,p_e3=zfse3u, cd_type='U', p_mask=umask, surf_crs=e2e3u, surf_msk_crs=e2e3u_msk) |
---|
483 | CALL crs_surf(p_e1=e1v, p_e2=e2v ,p_e3=zfse3v, cd_type='V', p_mask=vmask, surf_crs=e1e3v, surf_msk_crs=e1e3v_msk) |
---|
484 | |
---|
485 | |
---|
486 | !--------------------------------------------------------- |
---|
487 | ! 4. Coarse grid ocean volume and averaging weights |
---|
488 | !--------------------------------------------------------- |
---|
489 | ! 4.a. Ocean volume or area unmasked and masked |
---|
490 | |
---|
491 | !! ! jes. May not need ocean_volume_crs_t, ocean_volume_crs_w as calculated already in trc_init as cvol |
---|
492 | CALL crsfun( cd_type='T', cd_op='VOL', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_fse3=zfse3t, & |
---|
493 | & p_cfield3d_1=ocean_volume_crs_t, p_cfield3d_2=facvol_t ) |
---|
494 | |
---|
495 | r1_bt_crs(:,:,:) = 0._wp |
---|
496 | bt_crs(:,:,:) = ocean_volume_crs_t(:,:,:)* facvol_t(:,:,:) |
---|
497 | WHERE( bt_crs /= 0._wp ) r1_bt_crs(:,:,:) = 1._wp/bt_crs(:,:,:) |
---|
498 | |
---|
499 | CALL crsfun( cd_type='W', cd_op='VOL', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_fse3=zfse3w, & |
---|
500 | & p_cfield3d_1=ocean_volume_crs_w, p_cfield3d_2=facvol_w ) |
---|
501 | |
---|
502 | CALL crsfun( cd_type='U', cd_op='SUM', p_pmask=umask, p_e1=e1u, p_e2=e2u, p_fse3=zfse3u, p_cfield3d_1=facsurfu ) |
---|
503 | CALL crsfun( cd_type='V', cd_op='SUM', p_pmask=vmask, p_e1=e1v, p_e2=e2v, p_fse3=zfse3v, p_cfield3d_1=facsurfv ) |
---|
504 | |
---|
505 | |
---|
506 | ! 4.b. V, U and W surface area weights |
---|
507 | CALL crsfun( cd_type='U', cd_op='WGT', p_pmask=umask, p_e1=e1u, p_e2=e2u, p_fse3=zfse3u, p_cfield3d_1=crs_surfu_wgt ) |
---|
508 | CALL crsfun( cd_type='V', cd_op='WGT', p_pmask=vmask, p_e1=e1v, p_e2=e2v, p_fse3=zfse3v, p_cfield3d_1=crs_surfv_wgt ) |
---|
509 | CALL crsfun( cd_type='W', cd_op='WGT', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_cfield3d_1=crs_surfw_wgt ) |
---|
510 | |
---|
511 | ! 4.c. T volume weights |
---|
512 | CALL crsfun( cd_type='T', cd_op='WGT', p_pmask=tmask, p_e1=e1t, p_e2=e2t, p_fse3=zfse3t, p_cfield3d_1=crs_volt_wgt ) |
---|
513 | |
---|
514 | !--------------------------------------------------------- |
---|
515 | ! 5. Write out coarse meshmask (see OPA_SRC/DOM/domwri.F90 for ideas later) |
---|
516 | !--------------------------------------------------------- |
---|
517 | IF ( nn_msh_crs > 0 ) CALL crs_dom_wri |
---|
518 | |
---|
519 | WRITE(numout,*) 'crs_init done' |
---|
520 | |
---|
521 | !--------------------------------------------------------- |
---|
522 | ! 7. Finish and clean-up |
---|
523 | !--------------------------------------------------------- |
---|
524 | DEALLOCATE( zmbk ) |
---|
525 | DEALLOCATE( zfse3t, zfse3u, zfse3v, zfse3f ) |
---|
526 | DEALLOCATE( zfse3w, zfse3t_n, zfse3t_b ) |
---|
527 | |
---|
528 | |
---|
529 | END SUBROUTINE crs_init |
---|
530 | |
---|
531 | !!====================================================================== |
---|
532 | |
---|
533 | END MODULE crsini |
---|