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.

domain.F90 in branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 3820

Last change on this file since 3820 was 3764, checked in by smasson, 11 years ago
dev_MERGE_2012: report bugfixes done in the trunk from r3555 to r3763 into dev_MERGE_2012
Property svn:keywords set to `Id`
File size: 18.7 KB

Rev	Line
[3]	1	MODULE domain
	2	!!==============================================================================
	3	!! * MODULE domain *
	4	!! Ocean initialization : domain initialization
	5	!!==============================================================================
[1438]	6	!! History : OPA ! 1990-10 (C. Levy - G. Madec) Original code
	7	!! ! 1992-01 (M. Imbard) insert time step initialization
	8	!! ! 1996-06 (G. Madec) generalized vertical coordinate
	9	!! ! 1997-02 (G. Madec) creation of domwri.F
	10	!! ! 2001-05 (E.Durand - G. Madec) insert closed sea
	11	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
	12	!! 2.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization
[2528]	13	!! 3.3 ! 2010-11 (G. Madec) initialisation in C1D configuration
[3]	14	!!----------------------------------------------------------------------
[1438]	15
	16	!!----------------------------------------------------------------------
[3]	17	!! dom_init : initialize the space and time domain
	18	!! dom_nam : read and contral domain namelists
	19	!! dom_ctl : control print for the ocean domain
	20	!!----------------------------------------------------------------------
[2528]	21	USE oce ! ocean variables
	22	USE dom_oce ! domain: ocean
[888]	23	USE sbc_oce ! surface boundary condition: ocean
[719]	24	USE phycst ! physical constants
[1601]	25	USE closea ! closed seas
[719]	26	USE in_out_manager ! I/O manager
[3]	27	USE lib_mpp ! distributed memory computing library
	28
	29	USE domhgr ! domain: set the horizontal mesh
	30	USE domzgr ! domain: set the vertical mesh
	31	USE domstp ! domain: set the time-step
	32	USE dommsk ! domain: set the mask system
	33	USE domwri ! domain: write the meshmask file
[592]	34	USE domvvl ! variable volume
[2528]	35	USE c1d ! 1D vertical configuration
	36	USE dyncor_c1d ! Coriolis term (c1d case) (cor_c1d routine)
[3294]	37	USE timing ! Timing
[3680]	38	USE lbclnk ! ocean lateral boundary condition (or mpp link)
[3]	39
	40	IMPLICIT NONE
	41	PRIVATE
	42
[1438]	43	PUBLIC dom_init ! called by opa.F90
[3]	44
	45	!! * Substitutions
	46	# include "domzgr_substitute.h90"
[1438]	47	!!-------------------------------------------------------------------------
[2528]	48	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	49	!! $Id$
[2528]	50	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1438]	51	!!-------------------------------------------------------------------------
[3]	52	CONTAINS
	53
	54	SUBROUTINE dom_init
	55	!!----------------------------------------------------------------------
	56	!! * ROUTINE dom_init *
	57	!!
	58	!! ** Purpose : Domain initialization. Call the routines that are
[1601]	59	!! required to create the arrays which define the space
	60	!! and time domain of the ocean model.
[3]	61	!!
[1601]	62	!! ** Method : - dom_msk: compute the masks from the bathymetry file
	63	!! - dom_hgr: compute or read the horizontal grid-point position
	64	!! and scale factors, and the coriolis factor
	65	!! - dom_zgr: define the vertical coordinate and the bathymetry
	66	!! - dom_stp: defined the model time step
	67	!! - dom_wri: create the meshmask file if nmsh=1
[2528]	68	!! - 1D configuration, move Coriolis, u and v at T-point
[3]	69	!!----------------------------------------------------------------------
[3764]	70	INTEGER :: jk ! dummy loop argument
	71	INTEGER :: iconf = 0 ! local integers
[3]	72	!!----------------------------------------------------------------------
[1601]	73	!
[3764]	74	IF( nn_timing == 1 ) CALL timing_start('dom_init')
[3294]	75	!
[3]	76	IF(lwp) THEN
	77	WRITE(numout,*)
	78	WRITE(numout,*) 'dom_init : domain initialization'
	79	WRITE(numout,*) '~~~~~~~~'
	80	ENDIF
[1601]	81	!
	82	CALL dom_nam ! read namelist ( namrun, namdom, namcla )
	83	CALL dom_clo ! Closed seas and lake
	84	CALL dom_hgr ! Horizontal mesh
	85	CALL dom_zgr ! Vertical mesh and bathymetry
	86	CALL dom_msk ! Masks
[3680]	87	IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency
[1601]	88	IF( lk_vvl ) CALL dom_vvl ! Vertical variable mesh
	89	!
[3764]	90	IF( lk_c1d ) CALL cor_c1d ! 1D configuration: Coriolis set at T-point
[2528]	91	!
[3764]	92	hu(:,:) = 0._wp ! Ocean depth at U- and V-points
	93	hv(:,:) = 0._wp
[3]	94	DO jk = 1, jpk
	95	hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk)
	96	hv(:,:) = hv(:,:) + fse3v(:,:,jk) * vmask(:,:,jk)
	97	END DO
[1601]	98	! ! Inverse of the local depth
[3764]	99	hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)
	100	hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
[216]	101
[1601]	102	CALL dom_stp ! time step
	103	IF( nmsh /= 0 ) CALL dom_wri ! Create a domain file
	104	IF( .NOT.ln_rstart ) CALL dom_ctl ! Domain control
[1438]	105	!
[3764]	106	IF( nn_timing == 1 ) CALL timing_stop('dom_init')
[3294]	107	!
[3]	108	END SUBROUTINE dom_init
	109
	110
	111	SUBROUTINE dom_nam
	112	!!----------------------------------------------------------------------
	113	!! * ROUTINE dom_nam *
	114	!!
	115	!! ** Purpose : read domaine namelists and print the variables.
	116	!!
	117	!! ** input : - namrun namelist
	118	!! - namdom namelist
	119	!! - namcla namelist
[2528]	120	!! - namnc4 namelist ! "key_netcdf4" only
[3]	121	!!----------------------------------------------------------------------
	122	USE ioipsl
[1601]	123	NAMELIST/namrun/ nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl, &
	124	& nn_it000, nn_itend , nn_date0 , nn_leapy , nn_istate , nn_stock , &
	125	& nn_write, ln_dimgnnn, ln_mskland , ln_clobber , nn_chunksz
[2528]	126	NAMELIST/namdom/ nn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh , rn_hmin, &
	127	& nn_acc , rn_atfp , rn_rdt , rn_rdtmin , &
[1601]	128	& rn_rdtmax, rn_rdth , nn_baro , nn_closea
	129	NAMELIST/namcla/ nn_cla
[2528]	130	#if defined key_netcdf4
	131	NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
	132	#endif
[3]	133	!!----------------------------------------------------------------------
	134
[1601]	135	REWIND( numnam ) ! Namelist namrun : parameters of the run
	136	READ ( numnam, namrun )
	137	!
	138	IF(lwp) THEN ! control print
[3]	139	WRITE(numout,*)
	140	WRITE(numout,*) 'dom_nam : domain initialization through namelist read'
	141	WRITE(numout,*) '~~~~~~~ '
[1601]	142	WRITE(numout,*) ' Namelist namrun'
	143	WRITE(numout,*) ' job number nn_no = ', nn_no
	144	WRITE(numout,*) ' experiment name for output cn_exp = ', cn_exp
	145	WRITE(numout,*) ' restart logical ln_rstart = ', ln_rstart
[1604]	146	WRITE(numout,*) ' control of time step nn_rstctl = ', nn_rstctl
[1601]	147	WRITE(numout,*) ' number of the first time step nn_it000 = ', nn_it000
	148	WRITE(numout,*) ' number of the last time step nn_itend = ', nn_itend
	149	WRITE(numout,*) ' initial calendar date aammjj nn_date0 = ', nn_date0
	150	WRITE(numout,*) ' leap year calendar (0/1) nn_leapy = ', nn_leapy
	151	WRITE(numout,*) ' initial state output nn_istate = ', nn_istate
	152	WRITE(numout,*) ' frequency of restart file nn_stock = ', nn_stock
	153	WRITE(numout,*) ' frequency of output file nn_write = ', nn_write
	154	WRITE(numout,*) ' multi file dimgout ln_dimgnnn = ', ln_dimgnnn
	155	WRITE(numout,*) ' mask land points ln_mskland = ', ln_mskland
	156	WRITE(numout,*) ' overwrite an existing file ln_clobber = ', ln_clobber
	157	WRITE(numout,*) ' NetCDF chunksize (bytes) nn_chunksz = ', nn_chunksz
[3]	158	ENDIF
	159
[1601]	160	no = nn_no ! conversion DOCTOR names into model names (this should disappear soon)
	161	cexper = cn_exp
	162	nrstdt = nn_rstctl
	163	nit000 = nn_it000
	164	nitend = nn_itend
	165	ndate0 = nn_date0
	166	nleapy = nn_leapy
	167	ninist = nn_istate
	168	nstock = nn_stock
	169	nwrite = nn_write
[3]	170
[1601]	171
	172	! ! control of output frequency
[1335]	173	IF ( nstock == 0 .OR. nstock > nitend ) THEN
[1601]	174	WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
[783]	175	CALL ctl_warn( ctmp1 )
[1335]	176	nstock = nitend
[3]	177	ENDIF
	178	IF ( nwrite == 0 ) THEN
[1601]	179	WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
[783]	180	CALL ctl_warn( ctmp1 )
	181	nwrite = nitend
[3]	182	ENDIF
	183
[2528]	184	#if defined key_agrif
[1601]	185	IF( Agrif_Root() ) THEN
[2528]	186	#endif
	187	SELECT CASE ( nleapy ) ! Choose calendar for IOIPSL
	188	CASE ( 1 )
	189	CALL ioconf_calendar('gregorian')
	190	IF(lwp) WRITE(numout,*) ' The IOIPSL calendar is "gregorian", i.e. leap year'
	191	CASE ( 0 )
	192	CALL ioconf_calendar('noleap')
	193	IF(lwp) WRITE(numout,*) ' The IOIPSL calendar is "noleap", i.e. no leap year'
	194	CASE ( 30 )
	195	CALL ioconf_calendar('360d')
	196	IF(lwp) WRITE(numout,*) ' The IOIPSL calendar is "360d", i.e. 360 days in a year'
	197	END SELECT
	198	#if defined key_agrif
[1601]	199	ENDIF
[2528]	200	#endif
[3]	201
[2528]	202	REWIND( numnam ) ! Namelist namdom : space & time domain (bathymetry, mesh, timestep)
[3]	203	READ ( numnam, namdom )
	204
	205	IF(lwp) THEN
[72]	206	WRITE(numout,*)
[1601]	207	WRITE(numout,*) ' Namelist namdom : space & time domain'
	208	WRITE(numout,*) ' flag read/compute bathymetry nn_bathy = ', nn_bathy
[2528]	209	WRITE(numout,*) ' min depth of the ocean (>0) or rn_hmin = ', rn_hmin
	210	WRITE(numout,*) ' min number of ocean level (<0) '
[1601]	211	WRITE(numout,*) ' minimum thickness of partial rn_e3zps_min = ', rn_e3zps_min, ' (m)'
	212	WRITE(numout,*) ' step level rn_e3zps_rat = ', rn_e3zps_rat
	213	WRITE(numout,*) ' create mesh/mask file(s) nn_msh = ', nn_msh
[2528]	214	WRITE(numout,*) ' = 0 no file created '
	215	WRITE(numout,*) ' = 1 mesh_mask '
	216	WRITE(numout,*) ' = 2 mesh and mask '
	217	WRITE(numout,*) ' = 3 mesh_hgr, msh_zgr and mask'
[1601]	218	WRITE(numout,*) ' ocean time step rn_rdt = ', rn_rdt
	219	WRITE(numout,*) ' asselin time filter parameter rn_atfp = ', rn_atfp
	220	WRITE(numout,*) ' time-splitting: nb of sub time-step nn_baro = ', nn_baro
	221	WRITE(numout,*) ' acceleration of converge nn_acc = ', nn_acc
	222	WRITE(numout,*) ' nn_acc=1: surface tracer rdt rn_rdtmin = ', rn_rdtmin
	223	WRITE(numout,*) ' bottom tracer rdt rdtmax = ', rn_rdtmax
	224	WRITE(numout,*) ' depth of transition rn_rdth = ', rn_rdth
	225	WRITE(numout,*) ' suppression of closed seas (=0) nn_closea = ', nn_closea
[223]	226	ENDIF
	227
[1601]	228	ntopo = nn_bathy ! conversion DOCTOR names into model names (this should disappear soon)
	229	e3zps_min = rn_e3zps_min
	230	e3zps_rat = rn_e3zps_rat
	231	nmsh = nn_msh
	232	nacc = nn_acc
	233	atfp = rn_atfp
	234	rdt = rn_rdt
	235	rdtmin = rn_rdtmin
	236	rdtmax = rn_rdtmin
	237	rdth = rn_rdth
	238
[2528]	239	REWIND( numnam ) ! Namelist cross land advection
[3]	240	READ ( numnam, namcla )
	241	IF(lwp) THEN
[72]	242	WRITE(numout,*)
[1601]	243	WRITE(numout,*) ' Namelist namcla'
	244	WRITE(numout,*) ' cross land advection nn_cla = ', nn_cla
[3]	245	ENDIF
	246
[2528]	247	#if defined key_netcdf4
	248	! ! NetCDF 4 case ("key_netcdf4" defined)
	249	REWIND( numnam ) ! Namelist namnc4 : netcdf4 chunking parameters
	250	READ ( numnam, namnc4 )
	251	IF(lwp) THEN ! control print
	252	WRITE(numout,*)
	253	WRITE(numout,*) ' Namelist namnc4 - Netcdf4 chunking parameters'
	254	WRITE(numout,*) ' number of chunks in i-dimension nn_nchunks_i = ', nn_nchunks_i
	255	WRITE(numout,*) ' number of chunks in j-dimension nn_nchunks_j = ', nn_nchunks_j
	256	WRITE(numout,*) ' number of chunks in k-dimension nn_nchunks_k = ', nn_nchunks_k
	257	WRITE(numout,*) ' apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
	258	ENDIF
[1601]	259
[2528]	260	! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
	261	! Note the chunk size in the unlimited (time) dimension will be fixed at 1
	262	snc4set%ni = nn_nchunks_i
	263	snc4set%nj = nn_nchunks_j
	264	snc4set%nk = nn_nchunks_k
	265	snc4set%luse = ln_nc4zip
	266	#else
	267	snc4set%luse = .FALSE. ! No NetCDF 4 case
	268	#endif
[1438]	269	!
[3]	270	END SUBROUTINE dom_nam
	271
	272
	273	SUBROUTINE dom_ctl
	274	!!----------------------------------------------------------------------
	275	!! * ROUTINE dom_ctl *
	276	!!
	277	!! ** Purpose : Domain control.
	278	!!
	279	!! ** Method : compute and print extrema of masked scale factors
	280	!!----------------------------------------------------------------------
	281	INTEGER :: iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
[1601]	282	INTEGER, DIMENSION(2) :: iloc !
[3]	283	REAL(wp) :: ze1min, ze1max, ze2min, ze2max
	284	!!----------------------------------------------------------------------
[1601]	285	!
	286	IF(lk_mpp) THEN
[181]	287	CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 )
	288	CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 )
	289	CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 )
	290	CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 )
	291	ELSE
[3764]	292	ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
	293	ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
	294	ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
	295	ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
[32]	296
[3764]	297	iloc = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
[181]	298	iimi1 = iloc(1) + nimpp - 1
	299	ijmi1 = iloc(2) + njmpp - 1
[3764]	300	iloc = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
[181]	301	iimi2 = iloc(1) + nimpp - 1
	302	ijmi2 = iloc(2) + njmpp - 1
[3764]	303	iloc = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
[181]	304	iima1 = iloc(1) + nimpp - 1
	305	ijma1 = iloc(2) + njmpp - 1
[3764]	306	iloc = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
[181]	307	iima2 = iloc(1) + nimpp - 1
	308	ijma2 = iloc(2) + njmpp - 1
[32]	309	ENDIF
[3]	310	IF(lwp) THEN
[1601]	311	WRITE(numout,*)
	312	WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
	313	WRITE(numout,*) '~~~~~~~'
[181]	314	WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
	315	WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
	316	WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
	317	WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
[3]	318	ENDIF
[1438]	319	!
[3]	320	END SUBROUTINE dom_ctl
	321
[3680]	322	SUBROUTINE dom_stiff
	323	!!----------------------------------------------------------------------
	324	!! * ROUTINE dom_stiff *
	325	!!
	326	!! ** Purpose : Diagnose maximum grid stiffness/hydrostatic consistency
	327	!!
	328	!! ** Method : Compute Haney (1991) hydrostatic condition ratio
	329	!! Save the maximum in the vertical direction
	330	!! (this number is only relevant in s-coordinates)
	331	!!
	332	!! Haney, R. L., 1991: On the pressure gradient force
	333	!! over steep topography in sigma coordinate ocean models.
	334	!! J. Phys. Oceanogr., 21, 610???619.
	335	!!----------------------------------------------------------------------
	336	INTEGER :: ji, jj, jk
	337	REAL(wp) :: zrxmax
	338	REAL(wp), DIMENSION(4) :: zr1
	339	!!----------------------------------------------------------------------
	340	rx1(:,:) = 0.e0
	341	zrxmax = 0.e0
	342	zr1(:) = 0.e0
	343
	344	DO ji = 2, jpim1
	345	DO jj = 2, jpjm1
	346	DO jk = 1, jpkm1
	347	zr1(1) = umask(ji-1,jj ,jk) *abs( (gdepw(ji ,jj ,jk )-gdepw(ji-1,jj ,jk ) &
	348	& +gdepw(ji ,jj ,jk+1)-gdepw(ji-1,jj ,jk+1)) &
	349	& /(gdepw(ji ,jj ,jk )+gdepw(ji-1,jj ,jk ) &
	350	& -gdepw(ji ,jj ,jk+1)-gdepw(ji-1,jj ,jk+1) + rsmall) )
	351	zr1(2) = umask(ji ,jj ,jk) *abs( (gdepw(ji+1,jj ,jk )-gdepw(ji ,jj ,jk ) &
	352	& +gdepw(ji+1,jj ,jk+1)-gdepw(ji ,jj ,jk+1)) &
	353	& /(gdepw(ji+1,jj ,jk )+gdepw(ji ,jj ,jk ) &
	354	& -gdepw(ji+1,jj ,jk+1)-gdepw(ji ,jj ,jk+1) + rsmall) )
	355	zr1(3) = vmask(ji ,jj ,jk) *abs( (gdepw(ji ,jj+1,jk )-gdepw(ji ,jj ,jk ) &
	356	& +gdepw(ji ,jj+1,jk+1)-gdepw(ji ,jj ,jk+1)) &
	357	& /(gdepw(ji ,jj+1,jk )+gdepw(ji ,jj ,jk ) &
	358	& -gdepw(ji ,jj+1,jk+1)-gdepw(ji ,jj ,jk+1) + rsmall) )
	359	zr1(4) = vmask(ji ,jj-1,jk) *abs( (gdepw(ji ,jj ,jk )-gdepw(ji ,jj-1,jk ) &
	360	& +gdepw(ji ,jj ,jk+1)-gdepw(ji ,jj-1,jk+1)) &
	361	& /(gdepw(ji ,jj ,jk )+gdepw(ji ,jj-1,jk ) &
	362	& -gdepw(ji, jj ,jk+1)-gdepw(ji ,jj-1,jk+1) + rsmall) )
	363	zrxmax = MAXVAL(zr1(1:4))
	364	rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)
	365	END DO
	366	END DO
	367	END DO
	368
	369	CALL lbc_lnk( rx1, 'T', 1. )
	370
	371	zrxmax = MAXVAL(rx1)
	372
	373	IF( lk_mpp ) CALL mpp_max( zrxmax ) ! max over the global domain
	374
	375	IF(lwp) THEN
	376	WRITE(numout,*)
	377	WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
	378	WRITE(numout,*) '~~~~~~~~~'
	379	ENDIF
	380
	381	END SUBROUTINE dom_stiff
	382
	383
	384
[3]	385	!!======================================================================
	386	END MODULE domain

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

Download in other formats: