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

Context Navigation

source: branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsini.F90 @ 3779

Last change on this file since 3779 was 3779, checked in by cetlod, 11 years ago
2013/dev_r3411_CNRS4_IOCRS/NEMOGCM : lines are toot long for gfortran...
File size: 23.0 KB

Rev	Line
[3622]	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 CdA or CdV 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
[3735]	69	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zfse3t, zfse3u, zfse3v, zfse3f
	70	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zfse3w, zfse3t_n, zfse3t_b
[3622]	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
[3727]	96	WRITE(numout,*) ' nn_msh_crs = ', nn_msh_crs
[3622]	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-2jpreci + (jpni-1) ) / jpni + 2jpreci
	114	jpj_crs = ( jpjglo_crs-2jprecj + (jpnj-1) ) / jpnj + 2jprecj
	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
[3778]	195	!cc ijie = (ji*nn_factx)-nn_factx+1
	196	ijie = (ji*nn_factx)-nn_factx !cc
[3622]	197	ijis = ijie-nn_factx+1
	198
	199	IF ( ji == jpiglo_crsm1 ) THEN
[3778]	200	IF ( ((jpiglo-1)-ijie) <= nn_factx ) ijie = jpiglo-2 ! ijie = jpiglo-1 !cc
[3622]	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
[3778]	210	DO jj = 2, jpjglo_crsm1
	211	! cc ijje = ijjgloT-nn_facty*(jj-2)
	212	ijje = ijjgloT-nn_facty*(jj-2) - 1
[3622]	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
[3778]	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
[3622]	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
[3778]	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
[3622]	320	WRITE(numout,*) 'crsini. count 2'
	321
[3778]	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
[3622]	325
	326	WRITE(numout,*) 'crsini. count 3'
[3778]	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
[3622]	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
[3778]	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 )
[3622]	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
[3735]	441	ALLOCATE( zfse3t(jpi,jpj,jpk), zfse3u(jpi,jpj,jpk), zfse3v(jpi,jpj,jpk), zfse3f(jpi,jpj,jpk), &
[3622]	442	& zfse3w(jpi,jpj,jpk), zfse3t_n(jpi,jpj,jpk), zfse3t_b(jpi,jpj,jpk) )
	443	zfse3t(:,:,:) = fse3t(:,:,:)
	444	zfse3u(:,:,:) = fse3u(:,:,:)
	445	zfse3v(:,:,:) = fse3v(:,:,:)
[3735]	446	zfse3f(:,:,:) = fse3f(:,:,:)
[3622]	447	zfse3w(:,:,:) = fse3w(:,:,:)
[3778]	448
	449
[3622]	450
[3778]	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
[3622]	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)
[3735]	467	IF ( e3f_crs(ji,jj,jk) == 0 ) e3f_crs(ji,jj,jk) = e3t_0(jk)
[3622]	468	ENDDO
	469	ENDDO
	470	ENDDO
	471
	472	! 3.d.3 Vertical depth (meters)
	473
[3779]	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 )
[3622]	478
[3778]	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)
[3622]	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 )
[3778]	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(:,:,:)
[3622]	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	!---------------------------------------------------------
[3727]	517	IF ( nn_msh_crs > 0 ) CALL crs_dom_wri
[3622]	518
	519	WRITE(numout,*) 'crs_init done'
	520
	521	!---------------------------------------------------------
	522	! 7. Finish and clean-up
	523	!---------------------------------------------------------
	524	DEALLOCATE( zmbk )
[3735]	525	DEALLOCATE( zfse3t, zfse3u, zfse3v, zfse3f )
	526	DEALLOCATE( zfse3w, zfse3t_n, zfse3t_b )
[3622]	527
	528
	529	END SUBROUTINE crs_init
	530
	531	!!======================================================================
	532
	533	END MODULE crsini

Note: See TracBrowser for help on using the repository browser.

Download in other formats: