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/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/DOM – NEMO

Context Navigation

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/DOM/domwri.F90 @ 10380

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

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

Download in other formats: