Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

domwri.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 16.6 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
	20	USE wet_dry, ONLY : ln_wd, ht_wd
	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
[3294]	26	USE timing ! Timing
[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	!!----------------------------------------------------------------------
[7646]	37	!! NEMO/OPA 4.0 , NEMO Consortium (2016)
[1152]	38	!! $Id$
[2528]	39	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[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	!!
	52	!! ** Method : Write in a file all the arrays generated in routines
	53	!! domhgr, domzgr, and dommsk. Note: the file contain depends on
	54	!! the vertical coord. used (z-coord, partial steps, s-coord)
[7646]	55	!! MOD(nn_msh, 3) = 1 : 'mesh_mask.nc' file
[3]	56	!! = 2 : 'mesh.nc' and mask.nc' files
[1929]	57	!! = 0 : 'mesh_hgr.nc', 'mesh_zgr.nc' and
[3]	58	!! 'mask.nc' files
	59	!! For huge size domain, use option 2 or 3 depending on your
	60	!! vertical coordinate.
	61	!!
[7646]	62	!! if nn_msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
	63	!! if 3 < nn_msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays
[2528]	64	!! corresponding to the depth of the bottom t- and w-points
[7646]	65	!! if 6 < nn_msh <= 9: write 2D arrays corresponding to the depth and the
[2528]	66	!! thickness (e3[tw]_ps) of the bottom points
[1929]	67	!!
[2528]	68	!! ** output file : meshmask.nc : domain size, horizontal grid-point position,
	69	!! masks, depth and vertical scale factors
[3]	70	!!----------------------------------------------------------------------
[7646]	71	INTEGER :: inum ! temprary units for 'mesh_mask.nc' file
	72	CHARACTER(len=21) :: clnam ! filename (mesh and mask informations)
[2528]	73	INTEGER :: ji, jj, jk ! dummy loop indices
[7646]	74	INTEGER :: izco, izps, isco, icav
	75	!
[7910]	76	REAL(wp), DIMENSION(jpi,jpj) :: zprt, zprw ! 2D workspace
	77	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdepu, zdepv ! 3D workspace
[2715]	78	!!----------------------------------------------------------------------
[3294]	79	!
	80	IF( nn_timing == 1 ) CALL timing_start('dom_wri')
	81	!
	82	!
[1590]	83	IF(lwp) WRITE(numout,*)
	84	IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
	85	IF(lwp) WRITE(numout,*) '~~~~~~~'
	86
[7646]	87	clnam = 'mesh_mask' ! filename (mesh and mask informations)
[1590]	88
[7646]	89	! ! ============================
	90	! ! create 'mesh_mask.nc' file
	91	! ! ============================
	92	CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE., kiolib = jprstlib )
	93	!
	94	! ! global domain size
	95	CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 )
	96	CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 )
	97	CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpkglo, wp), ktype = jp_i4 )
	98
	99	! ! domain characteristics
	100	CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
	101	! ! type of vertical coordinate
	102	IF( ln_zco ) THEN ; izco = 1 ; ELSE ; izco = 0 ; ENDIF
	103	IF( ln_zps ) THEN ; izps = 1 ; ELSE ; izps = 0 ; ENDIF
	104	IF( ln_sco ) THEN ; isco = 1 ; ELSE ; isco = 0 ; ENDIF
	105	CALL iom_rstput( 0, 0, inum, 'ln_zco' , REAL( izco, wp), ktype = jp_i4 )
	106	CALL iom_rstput( 0, 0, inum, 'ln_zps' , REAL( izps, wp), ktype = jp_i4 )
	107	CALL iom_rstput( 0, 0, inum, 'ln_sco' , REAL( isco, wp), ktype = jp_i4 )
	108	! ! ocean cavities under iceshelves
	109	IF( ln_isfcav ) THEN ; icav = 1 ; ELSE ; icav = 0 ; ENDIF
	110	CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
	111
	112	! ! masks
	113	CALL iom_rstput( 0, 0, inum, 'tmask', tmask, ktype = jp_i1 ) ! ! land-sea mask
	114	CALL iom_rstput( 0, 0, inum, 'umask', umask, ktype = jp_i1 )
	115	CALL iom_rstput( 0, 0, inum, 'vmask', vmask, ktype = jp_i1 )
	116	CALL iom_rstput( 0, 0, inum, 'fmask', fmask, ktype = jp_i1 )
[1161]	117
[2715]	118	CALL dom_uniq( zprw, 'T' )
[4990]	119	DO jj = 1, jpj
	120	DO ji = 1, jpi
[7646]	121	zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj) ! ! unique point mask
[4990]	122	END DO
	123	END DO ! ! unique point mask
[7646]	124	CALL iom_rstput( 0, 0, inum, 'tmaskutil', zprt, ktype = jp_i1 )
[2715]	125	CALL dom_uniq( zprw, 'U' )
[4990]	126	DO jj = 1, jpj
	127	DO ji = 1, jpi
[7646]	128	zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj) ! ! unique point mask
[4990]	129	END DO
	130	END DO
[7646]	131	CALL iom_rstput( 0, 0, inum, 'umaskutil', zprt, ktype = jp_i1 )
[2715]	132	CALL dom_uniq( zprw, 'V' )
[4990]	133	DO jj = 1, jpj
	134	DO ji = 1, jpi
[7646]	135	zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj) ! ! unique point mask
[4990]	136	END DO
	137	END DO
[7646]	138	CALL iom_rstput( 0, 0, inum, 'vmaskutil', zprt, ktype = jp_i1 )
	139	!!gm ssfmask has been removed ==>> find another solution to defined fmaskutil
	140	!! Here we just remove the output of fmaskutil.
	141	! CALL dom_uniq( zprw, 'F' )
	142	! DO jj = 1, jpj
	143	! DO ji = 1, jpi
	144	! zprt(ji,jj) = ssfmask(ji,jj) * zprw(ji,jj) ! ! unique point mask
	145	! END DO
	146	! END DO
	147	! CALL iom_rstput( 0, 0, inum, 'fmaskutil', zprt, ktype = jp_i1 )
	148	!!gm
[1161]	149
	150	! ! horizontal mesh (inum3)
[7646]	151	CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 ) ! ! latitude
	152	CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 )
	153	CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 )
	154	CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 )
[1161]	155
[7646]	156	CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 ) ! ! longitude
	157	CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 )
	158	CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 )
	159	CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 )
[1161]	160
[7646]	161	CALL iom_rstput( 0, 0, inum, 'e1t', e1t, ktype = jp_r8 ) ! ! e1 scale factors
	162	CALL iom_rstput( 0, 0, inum, 'e1u', e1u, ktype = jp_r8 )
	163	CALL iom_rstput( 0, 0, inum, 'e1v', e1v, ktype = jp_r8 )
	164	CALL iom_rstput( 0, 0, inum, 'e1f', e1f, ktype = jp_r8 )
[1161]	165
[7646]	166	CALL iom_rstput( 0, 0, inum, 'e2t', e2t, ktype = jp_r8 ) ! ! e2 scale factors
	167	CALL iom_rstput( 0, 0, inum, 'e2u', e2u, ktype = jp_r8 )
	168	CALL iom_rstput( 0, 0, inum, 'e2v', e2v, ktype = jp_r8 )
	169	CALL iom_rstput( 0, 0, inum, 'e2f', e2f, ktype = jp_r8 )
[1161]	170
[7646]	171	CALL iom_rstput( 0, 0, inum, 'ff_f', ff_f, ktype = jp_r8 ) ! ! coriolis factor
	172	CALL iom_rstput( 0, 0, inum, 'ff_t', ff_t, ktype = jp_r8 )
[1161]	173
[2528]	174	! note that mbkt is set to 1 over land ==> use surface tmask
[4990]	175	zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp )
[7646]	176	CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 ) ! ! nb of ocean T-points
[4990]	177	zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
[7646]	178	CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 ) ! ! nb of ocean T-points
[4990]	179	zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
[7646]	180	CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 ) ! ! nb of ocean T-points
	181	! ! vertical mesh
	182	CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8 ) ! ! scale factors
	183	CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0, ktype = jp_r8 )
	184	CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0, ktype = jp_r8 )
	185	CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0, ktype = jp_r8 )
	186	!
	187	CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 ) ! stretched system
	188	CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d , ktype = jp_r8 )
	189	CALL iom_rstput( 0, 0, inum, 'gdept_0' , gdept_0 , ktype = jp_r8 )
	190	CALL iom_rstput( 0, 0, inum, 'gdepw_0' , gdepw_0 , ktype = jp_r8 )
	191	!
	192	IF( ln_sco ) THEN ! s-coordinate stiffness
	193	CALL dom_stiff( zprt )
	194	CALL iom_rstput( 0, 0, inum, 'stiffness', zprt ) ! Max. grid stiffness ratio
[1161]	195	ENDIF
[7646]	196	!
	197	IF( ln_wd ) THEN ! wetting and drying domain
	198	CALL iom_rstput( 0, 0, inum, 'ht_0' , ht_0 , ktype = jp_r8 )
	199	CALL iom_rstput( 0, 0, inum, 'ht_wd' , ht_wd , ktype = jp_r8 )
[1161]	200	ENDIF
	201	! ! ============================
[7646]	202	CALL iom_close( inum ) ! close the files
[1161]	203	! ! ============================
[2528]	204	!
[2715]	205	!
[3294]	206	IF( nn_timing == 1 ) CALL timing_stop('dom_wri')
	207	!
[1161]	208	END SUBROUTINE dom_wri
[3]	209
	210
[2715]	211	SUBROUTINE dom_uniq( puniq, cdgrd )
[1161]	212	!!----------------------------------------------------------------------
	213	!! * ROUTINE dom_uniq *
	214	!!
	215	!! ** Purpose : identify unique point of a grid (TUVF)
	216	!!
	217	!! ** Method : 1) aplly lbc_lnk on an array with different values for each element
	218	!! 2) check which elements have been changed
	219	!!----------------------------------------------------------------------
[2715]	220	CHARACTER(len=1) , INTENT(in ) :: cdgrd !
	221	REAL(wp), DIMENSION(:,:), INTENT(inout) :: puniq !
	222	!
	223	REAL(wp) :: zshift ! shift value link to the process number
	224	INTEGER :: ji ! dummy loop indices
	225	LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) :: lldbl ! store whether each point is unique or not
[7910]	226	REAL(wp), DIMENSION(jpi,jpj) :: ztstref
[1161]	227	!!----------------------------------------------------------------------
[3294]	228	!
	229	IF( nn_timing == 1 ) CALL timing_start('dom_uniq')
	230	!
	231	!
[1161]	232	! build an array with different values for each element
	233	! in mpp: make sure that these values are different even between process
	234	! -> apply a shift value according to the process number
	235	zshift = jpi * jpj * ( narea - 1 )
[2528]	236	ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
	237	!
[1161]	238	puniq(:,:) = ztstref(:,:) ! default definition
	239	CALL lbc_lnk( puniq, cdgrd, 1. ) ! apply boundary conditions
	240	lldbl(:,:,1) = puniq(:,:) == ztstref(:,:) ! check which values have been changed
[2528]	241	!
[1161]	242	puniq(:,:) = 1. ! default definition
	243	! fill only the inner part of the cpu with llbl converted into real
[2528]	244	puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
	245	!
[2715]	246	!
[3294]	247	IF( nn_timing == 1 ) CALL timing_stop('dom_uniq')
	248	!
[2715]	249	END SUBROUTINE dom_uniq
[3]	250
[7646]	251
	252	SUBROUTINE dom_stiff( px1 )
	253	!!----------------------------------------------------------------------
	254	!! * ROUTINE dom_stiff *
	255	!!
	256	!! ** Purpose : Diagnose maximum grid stiffness/hydrostatic consistency
	257	!!
	258	!! ** Method : Compute Haney (1991) hydrostatic condition ratio
	259	!! Save the maximum in the vertical direction
	260	!! (this number is only relevant in s-coordinates)
	261	!!
	262	!! Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
	263	!!----------------------------------------------------------------------
	264	REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL :: px1 ! stiffness
	265	!
	266	INTEGER :: ji, jj, jk
	267	REAL(wp) :: zrxmax
	268	REAL(wp), DIMENSION(4) :: zr1
	269	REAL(wp), DIMENSION(jpi,jpj) :: zx1
	270	!!----------------------------------------------------------------------
	271	zx1(:,:) = 0._wp
	272	zrxmax = 0._wp
	273	zr1(:) = 0._wp
	274	!
	275	DO ji = 2, jpim1
	276	DO jj = 2, jpjm1
	277	DO jk = 1, jpkm1
	278	!!gm remark: dk(gdepw) = e3t ===>>> possible simplification of the following calculation....
	279	!! especially since it is gde3w which is used to compute the pressure gradient
	280	!! furthermore, I think gdept_0 should be used below instead of w point in the numerator
	281	!! so that the ratio is computed at the same point (i.e. uw and vw) ....
	282	zr1(1) = ABS( ( gdepw_0(ji ,jj,jk )-gdepw_0(ji-1,jj,jk ) &
	283	& +gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) ) &
	284	& / ( gdepw_0(ji ,jj,jk )+gdepw_0(ji-1,jj,jk ) &
	285	& -gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall ) ) * umask(ji-1,jj,jk)
	286	zr1(2) = ABS( ( gdepw_0(ji+1,jj,jk )-gdepw_0(ji ,jj,jk ) &
	287	& +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) ) &
	288	& / ( gdepw_0(ji+1,jj,jk )+gdepw_0(ji ,jj,jk ) &
	289	& -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) + rsmall ) ) * umask(ji ,jj,jk)
	290	zr1(3) = ABS( ( gdepw_0(ji,jj+1,jk )-gdepw_0(ji,jj ,jk ) &
	291	& +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) ) &
	292	& / ( gdepw_0(ji,jj+1,jk )+gdepw_0(ji,jj ,jk ) &
	293	& -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) + rsmall ) ) * vmask(ji,jj ,jk)
	294	zr1(4) = ABS( ( gdepw_0(ji,jj ,jk )-gdepw_0(ji,jj-1,jk ) &
	295	& +gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) ) &
	296	& / ( gdepw_0(ji,jj ,jk )+gdepw_0(ji,jj-1,jk ) &
	297	& -gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall ) ) * vmask(ji,jj-1,jk)
	298	zrxmax = MAXVAL( zr1(1:4) )
	299	zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
	300	END DO
	301	END DO
	302	END DO
	303	CALL lbc_lnk( zx1, 'T', 1. )
	304	!
	305	IF( PRESENT( px1 ) ) px1 = zx1
	306	!
	307	zrxmax = MAXVAL( zx1 )
	308	!
	309	IF( lk_mpp ) CALL mpp_max( zrxmax ) ! max over the global domain
	310	!
	311	IF(lwp) THEN
	312	WRITE(numout,*)
	313	WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
	314	WRITE(numout,*) '~~~~~~~~~'
	315	ENDIF
	316	!
	317	END SUBROUTINE dom_stiff
	318
[3]	319	!!======================================================================
	320	END MODULE domwri

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 7910

Download in other formats: