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.

domwri.F90 in NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM – NEMO

Context Navigation

source: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domwri.F90 @ 12143

Last change on this file since 12143 was 12143, checked in by mathiot, 4 years ago
update ENHANCE-02_ISF_nemo to 12072 (sette in progress)
Property svn:keywords set to `Id`
File size: 15.4 KB

Rev	Line
[3]	1	MODULE domwri
	2	!!======================================================================
	3	!! * MODULE domwri *
[2715]	4	!! Ocean initialization : write the ocean domain mesh file(s)
[3]	5	!!======================================================================
[2528]	6	!! History : OPA ! 1997-02 (G. Madec) Original code
	7	!! 8.1 ! 1999-11 (M. Imbard) NetCDF FORMAT with IOIPSL
	8	!! NEMO 1.0 ! 2002-08 (G. Madec) F90 and several file
[5836]	9	!! 3.0 ! 2008-01 (S. Masson) add dom_uniq
[7646]	10	!! 4.0 ! 2016-01 (G. Madec) simplified mesh_mask.nc file
[2528]	11	!!----------------------------------------------------------------------
[3]	12
	13	!!----------------------------------------------------------------------
[84]	14	!! dom_wri : create and write mesh and mask file(s)
[5836]	15	!! dom_uniq : identify unique point of a grid (TUVF)
[7646]	16	!! dom_stiff : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
[3]	17	!!----------------------------------------------------------------------
[11987]	18	!
[3]	19	USE dom_oce ! ocean space and time domain
[7646]	20	USE phycst , ONLY : rsmall
[9023]	21	USE wet_dry, ONLY : ll_wd ! Wetting and drying
[7646]	22	!
[2528]	23	USE in_out_manager ! I/O manager
	24	USE iom ! I/O library
	25	USE lbclnk ! lateral boundary conditions - mpp exchanges
	26	USE lib_mpp ! MPP library
[3]	27
	28	IMPLICIT NONE
	29	PRIVATE
	30
[5836]	31	PUBLIC dom_wri ! routine called by inidom.F90
[7646]	32	PUBLIC dom_stiff ! routine called by inidom.F90
	33
[1590]	34	!! * Substitutions
	35	# include "vectopt_loop_substitute.h90"
[3]	36	!!----------------------------------------------------------------------
[9598]	37	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]	38	!! $Id$
[10068]	39	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	40	!!----------------------------------------------------------------------
	41	CONTAINS
	42
	43	SUBROUTINE dom_wri
	44	!!----------------------------------------------------------------------
	45	!! * ROUTINE dom_wri *
	46	!!
	47	!! ** Purpose : Create the NetCDF file(s) which contain(s) all the
	48	!! ocean domain informations (mesh and mask arrays). This (these)
	49	!! file(s) is (are) used for visualisation (SAXO software) and
	50	!! diagnostic computation.
	51	!!
[9169]	52	!! ** Method : create a file with all domain related arrays
[3]	53	!!
[2528]	54	!! ** output file : meshmask.nc : domain size, horizontal grid-point position,
	55	!! masks, depth and vertical scale factors
[3]	56	!!----------------------------------------------------------------------
[7646]	57	INTEGER :: inum ! temprary units for 'mesh_mask.nc' file
	58	CHARACTER(len=21) :: clnam ! filename (mesh and mask informations)
[2528]	59	INTEGER :: ji, jj, jk ! dummy loop indices
[7646]	60	INTEGER :: izco, izps, isco, icav
	61	!
[9019]	62	REAL(wp), DIMENSION(jpi,jpj) :: zprt, zprw ! 2D workspace
	63	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdepu, zdepv ! 3D workspace
[2715]	64	!!----------------------------------------------------------------------
[3294]	65	!
[1590]	66	IF(lwp) WRITE(numout,*)
	67	IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
	68	IF(lwp) WRITE(numout,*) '~~~~~~~'
	69
[7646]	70	clnam = 'mesh_mask' ! filename (mesh and mask informations)
[1590]	71
[7646]	72	! ! ============================
	73	! ! create 'mesh_mask.nc' file
	74	! ! ============================
[10425]	75	CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE. )
[7646]	76	!
	77	! ! global domain size
	78	CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 )
	79	CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 )
	80	CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpkglo, wp), ktype = jp_i4 )
	81
	82	! ! domain characteristics
	83	CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
	84	! ! type of vertical coordinate
	85	IF( ln_zco ) THEN ; izco = 1 ; ELSE ; izco = 0 ; ENDIF
	86	IF( ln_zps ) THEN ; izps = 1 ; ELSE ; izps = 0 ; ENDIF
	87	IF( ln_sco ) THEN ; isco = 1 ; ELSE ; isco = 0 ; ENDIF
	88	CALL iom_rstput( 0, 0, inum, 'ln_zco' , REAL( izco, wp), ktype = jp_i4 )
	89	CALL iom_rstput( 0, 0, inum, 'ln_zps' , REAL( izps, wp), ktype = jp_i4 )
	90	CALL iom_rstput( 0, 0, inum, 'ln_sco' , REAL( isco, wp), ktype = jp_i4 )
	91	! ! ocean cavities under iceshelves
	92	IF( ln_isfcav ) THEN ; icav = 1 ; ELSE ; icav = 0 ; ENDIF
	93	CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
	94
	95	! ! masks
	96	CALL iom_rstput( 0, 0, inum, 'tmask', tmask, ktype = jp_i1 ) ! ! land-sea mask
	97	CALL iom_rstput( 0, 0, inum, 'umask', umask, ktype = jp_i1 )
	98	CALL iom_rstput( 0, 0, inum, 'vmask', vmask, ktype = jp_i1 )
	99	CALL iom_rstput( 0, 0, inum, 'fmask', fmask, ktype = jp_i1 )
[1161]	100
[2715]	101	CALL dom_uniq( zprw, 'T' )
[4990]	102	DO jj = 1, jpj
	103	DO ji = 1, jpi
[7646]	104	zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj) ! ! unique point mask
[4990]	105	END DO
	106	END DO ! ! unique point mask
[7646]	107	CALL iom_rstput( 0, 0, inum, 'tmaskutil', zprt, ktype = jp_i1 )
[2715]	108	CALL dom_uniq( zprw, 'U' )
[4990]	109	DO jj = 1, jpj
	110	DO ji = 1, jpi
[7646]	111	zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj) ! ! unique point mask
[4990]	112	END DO
	113	END DO
[7646]	114	CALL iom_rstput( 0, 0, inum, 'umaskutil', zprt, ktype = jp_i1 )
[2715]	115	CALL dom_uniq( zprw, 'V' )
[4990]	116	DO jj = 1, jpj
	117	DO ji = 1, jpi
[7646]	118	zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj) ! ! unique point mask
[4990]	119	END DO
	120	END DO
[7646]	121	CALL iom_rstput( 0, 0, inum, 'vmaskutil', zprt, ktype = jp_i1 )
	122	!!gm ssfmask has been removed ==>> find another solution to defined fmaskutil
	123	!! Here we just remove the output of fmaskutil.
	124	! CALL dom_uniq( zprw, 'F' )
	125	! DO jj = 1, jpj
	126	! DO ji = 1, jpi
	127	! zprt(ji,jj) = ssfmask(ji,jj) * zprw(ji,jj) ! ! unique point mask
	128	! END DO
	129	! END DO
	130	! CALL iom_rstput( 0, 0, inum, 'fmaskutil', zprt, ktype = jp_i1 )
	131	!!gm
[1161]	132
	133	! ! horizontal mesh (inum3)
[7646]	134	CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 ) ! ! latitude
	135	CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 )
	136	CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 )
	137	CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 )
[1161]	138
[7646]	139	CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 ) ! ! longitude
	140	CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 )
	141	CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 )
	142	CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 )
[1161]	143
[7646]	144	CALL iom_rstput( 0, 0, inum, 'e1t', e1t, ktype = jp_r8 ) ! ! e1 scale factors
	145	CALL iom_rstput( 0, 0, inum, 'e1u', e1u, ktype = jp_r8 )
	146	CALL iom_rstput( 0, 0, inum, 'e1v', e1v, ktype = jp_r8 )
	147	CALL iom_rstput( 0, 0, inum, 'e1f', e1f, ktype = jp_r8 )
[1161]	148
[7646]	149	CALL iom_rstput( 0, 0, inum, 'e2t', e2t, ktype = jp_r8 ) ! ! e2 scale factors
	150	CALL iom_rstput( 0, 0, inum, 'e2u', e2u, ktype = jp_r8 )
	151	CALL iom_rstput( 0, 0, inum, 'e2v', e2v, ktype = jp_r8 )
	152	CALL iom_rstput( 0, 0, inum, 'e2f', e2f, ktype = jp_r8 )
[1161]	153
[7646]	154	CALL iom_rstput( 0, 0, inum, 'ff_f', ff_f, ktype = jp_r8 ) ! ! coriolis factor
	155	CALL iom_rstput( 0, 0, inum, 'ff_t', ff_t, ktype = jp_r8 )
[1161]	156
[2528]	157	! note that mbkt is set to 1 over land ==> use surface tmask
[11395]	158	zprt(:,:) = REAL( mbkt(:,:) , wp )
[7646]	159	CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 ) ! ! nb of ocean T-points
[11395]	160	zprt(:,:) = REAL( mikt(:,:) , wp )
[7646]	161	CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 ) ! ! nb of ocean T-points
	162	! ! vertical mesh
[12143]	163	CALL iom_rstput( 0, 0, inum, 'e3t_1d', e3t_1d, ktype = jp_r8 ) ! ! scale factors
	164	CALL iom_rstput( 0, 0, inum, 'e3w_1d', e3w_1d, ktype = jp_r8 )
	165
	166	CALL iom_rstput( 0, 0, inum, 'e3t_0' , e3t_0 , ktype = jp_r8 )
	167	CALL iom_rstput( 0, 0, inum, 'e3u_0' , e3u_0 , ktype = jp_r8 )
	168	CALL iom_rstput( 0, 0, inum, 'e3v_0' , e3v_0 , ktype = jp_r8 )
	169	CALL iom_rstput( 0, 0, inum, 'e3f_0' , e3f_0 , ktype = jp_r8 )
	170	CALL iom_rstput( 0, 0, inum, 'e3w_0' , e3w_0 , ktype = jp_r8 )
	171	CALL iom_rstput( 0, 0, inum, 'e3uw_0', e3uw_0, ktype = jp_r8 )
	172	CALL iom_rstput( 0, 0, inum, 'e3vw_0', e3vw_0, ktype = jp_r8 )
[7646]	173	!
	174	CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 ) ! stretched system
	175	CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d , ktype = jp_r8 )
	176	CALL iom_rstput( 0, 0, inum, 'gdept_0' , gdept_0 , ktype = jp_r8 )
	177	CALL iom_rstput( 0, 0, inum, 'gdepw_0' , gdepw_0 , ktype = jp_r8 )
	178	!
	179	IF( ln_sco ) THEN ! s-coordinate stiffness
	180	CALL dom_stiff( zprt )
	181	CALL iom_rstput( 0, 0, inum, 'stiffness', zprt ) ! Max. grid stiffness ratio
[1161]	182	ENDIF
[7646]	183	!
[9023]	184	IF( ll_wd ) CALL iom_rstput( 0, 0, inum, 'ht_0' , ht_0 , ktype = jp_r8 )
	185
[1161]	186	! ! ============================
[7646]	187	CALL iom_close( inum ) ! close the files
[1161]	188	! ! ============================
	189	END SUBROUTINE dom_wri
[3]	190
	191
[2715]	192	SUBROUTINE dom_uniq( puniq, cdgrd )
[1161]	193	!!----------------------------------------------------------------------
	194	!! * ROUTINE dom_uniq *
	195	!!
	196	!! ** Purpose : identify unique point of a grid (TUVF)
	197	!!
	198	!! ** Method : 1) aplly lbc_lnk on an array with different values for each element
	199	!! 2) check which elements have been changed
	200	!!----------------------------------------------------------------------
[2715]	201	CHARACTER(len=1) , INTENT(in ) :: cdgrd !
	202	REAL(wp), DIMENSION(:,:), INTENT(inout) :: puniq !
	203	!
	204	REAL(wp) :: zshift ! shift value link to the process number
	205	INTEGER :: ji ! dummy loop indices
	206	LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) :: lldbl ! store whether each point is unique or not
[9019]	207	REAL(wp), DIMENSION(jpi,jpj) :: ztstref
[1161]	208	!!----------------------------------------------------------------------
[3294]	209	!
[1161]	210	! build an array with different values for each element
	211	! in mpp: make sure that these values are different even between process
	212	! -> apply a shift value according to the process number
	213	zshift = jpi * jpj * ( narea - 1 )
[2528]	214	ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
	215	!
[1161]	216	puniq(:,:) = ztstref(:,:) ! default definition
[10425]	217	CALL lbc_lnk( 'domwri', puniq, cdgrd, 1. ) ! apply boundary conditions
[1161]	218	lldbl(:,:,1) = puniq(:,:) == ztstref(:,:) ! check which values have been changed
[2528]	219	!
[1161]	220	puniq(:,:) = 1. ! default definition
	221	! fill only the inner part of the cpu with llbl converted into real
[2528]	222	puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
	223	!
[2715]	224	END SUBROUTINE dom_uniq
[3]	225
[7646]	226
	227	SUBROUTINE dom_stiff( px1 )
	228	!!----------------------------------------------------------------------
	229	!! * ROUTINE dom_stiff *
	230	!!
	231	!! ** Purpose : Diagnose maximum grid stiffness/hydrostatic consistency
	232	!!
	233	!! ** Method : Compute Haney (1991) hydrostatic condition ratio
	234	!! Save the maximum in the vertical direction
	235	!! (this number is only relevant in s-coordinates)
	236	!!
	237	!! Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
	238	!!----------------------------------------------------------------------
	239	REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL :: px1 ! stiffness
	240	!
	241	INTEGER :: ji, jj, jk
	242	REAL(wp) :: zrxmax
	243	REAL(wp), DIMENSION(4) :: zr1
	244	REAL(wp), DIMENSION(jpi,jpj) :: zx1
	245	!!----------------------------------------------------------------------
	246	zx1(:,:) = 0._wp
	247	zrxmax = 0._wp
	248	zr1(:) = 0._wp
	249	!
	250	DO ji = 2, jpim1
	251	DO jj = 2, jpjm1
	252	DO jk = 1, jpkm1
	253	!!gm remark: dk(gdepw) = e3t ===>>> possible simplification of the following calculation....
	254	!! especially since it is gde3w which is used to compute the pressure gradient
	255	!! furthermore, I think gdept_0 should be used below instead of w point in the numerator
	256	!! so that the ratio is computed at the same point (i.e. uw and vw) ....
	257	zr1(1) = ABS( ( gdepw_0(ji ,jj,jk )-gdepw_0(ji-1,jj,jk ) &
	258	& +gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) ) &
	259	& / ( gdepw_0(ji ,jj,jk )+gdepw_0(ji-1,jj,jk ) &
	260	& -gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall ) ) * umask(ji-1,jj,jk)
	261	zr1(2) = ABS( ( gdepw_0(ji+1,jj,jk )-gdepw_0(ji ,jj,jk ) &
	262	& +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) ) &
	263	& / ( gdepw_0(ji+1,jj,jk )+gdepw_0(ji ,jj,jk ) &
	264	& -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) + rsmall ) ) * umask(ji ,jj,jk)
	265	zr1(3) = ABS( ( gdepw_0(ji,jj+1,jk )-gdepw_0(ji,jj ,jk ) &
	266	& +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) ) &
	267	& / ( gdepw_0(ji,jj+1,jk )+gdepw_0(ji,jj ,jk ) &
	268	& -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) + rsmall ) ) * vmask(ji,jj ,jk)
	269	zr1(4) = ABS( ( gdepw_0(ji,jj ,jk )-gdepw_0(ji,jj-1,jk ) &
	270	& +gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) ) &
	271	& / ( gdepw_0(ji,jj ,jk )+gdepw_0(ji,jj-1,jk ) &
	272	& -gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall ) ) * vmask(ji,jj-1,jk)
	273	zrxmax = MAXVAL( zr1(1:4) )
	274	zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
	275	END DO
	276	END DO
	277	END DO
[10425]	278	CALL lbc_lnk( 'domwri', zx1, 'T', 1. )
[7646]	279	!
	280	IF( PRESENT( px1 ) ) px1 = zx1
	281	!
	282	zrxmax = MAXVAL( zx1 )
	283	!
[10425]	284	CALL mpp_max( 'domwri', zrxmax ) ! max over the global domain
[7646]	285	!
	286	IF(lwp) THEN
	287	WRITE(numout,*)
	288	WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
	289	WRITE(numout,*) '~~~~~~~~~'
	290	ENDIF
	291	!
	292	END SUBROUTINE dom_stiff
	293
[3]	294	!!======================================================================
	295	END MODULE domwri

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

Download in other formats: