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 branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 4306

Last change on this file since 4306 was 4292, checked in by cetlod, 11 years ago
dev_MERGE_2013 : 1st step of the merge, see ticket #1185
Property svn:keywords set to `Id`
File size: 16.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
[2715]	9	!! 3.0 ! 2008-01 (S. Masson) add dom_uniq
	10	!! 4.0 ! 2011-01 (A. R. Porter, STFC Daresbury) dynamical allocation
[2528]	11	!!----------------------------------------------------------------------
[3]	12
	13	!!----------------------------------------------------------------------
[84]	14	!! dom_wri : create and write mesh and mask file(s)
[2715]	15	!! dom_uniq :
[3]	16	!!----------------------------------------------------------------------
	17	USE dom_oce ! ocean space and time domain
[2528]	18	USE in_out_manager ! I/O manager
	19	USE iom ! I/O library
	20	USE lbclnk ! lateral boundary conditions - mpp exchanges
	21	USE lib_mpp ! MPP library
[3294]	22	USE wrk_nemo ! Memory allocation
	23	USE timing ! Timing
[3]	24
	25	IMPLICIT NONE
	26	PRIVATE
	27
	28	PUBLIC dom_wri ! routine called by inidom.F90
[1590]	29
	30	!! * Substitutions
	31	# include "vectopt_loop_substitute.h90"
[3]	32	!!----------------------------------------------------------------------
[2528]	33	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1152]	34	!! $Id$
[2528]	35	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	36	!!----------------------------------------------------------------------
	37	CONTAINS
	38
	39	SUBROUTINE dom_wri
	40	!!----------------------------------------------------------------------
	41	!! * ROUTINE dom_wri *
	42	!!
	43	!! ** Purpose : Create the NetCDF file(s) which contain(s) all the
	44	!! ocean domain informations (mesh and mask arrays). This (these)
	45	!! file(s) is (are) used for visualisation (SAXO software) and
	46	!! diagnostic computation.
	47	!!
	48	!! ** Method : Write in a file all the arrays generated in routines
	49	!! domhgr, domzgr, and dommsk. Note: the file contain depends on
	50	!! the vertical coord. used (z-coord, partial steps, s-coord)
[1929]	51	!! MOD(nmsh, 3) = 1 : 'mesh_mask.nc' file
[3]	52	!! = 2 : 'mesh.nc' and mask.nc' files
[1929]	53	!! = 0 : 'mesh_hgr.nc', 'mesh_zgr.nc' and
[3]	54	!! 'mask.nc' files
	55	!! For huge size domain, use option 2 or 3 depending on your
	56	!! vertical coordinate.
	57	!!
[1929]	58	!! if nmsh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
	59	!! if 3 < nmsh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays
[2528]	60	!! corresponding to the depth of the bottom t- and w-points
[1929]	61	!! if 6 < nmsh <= 9: write 2D arrays corresponding to the depth and the
[2528]	62	!! thickness (e3[tw]_ps) of the bottom points
[1929]	63	!!
[2528]	64	!! ** output file : meshmask.nc : domain size, horizontal grid-point position,
	65	!! masks, depth and vertical scale factors
[3]	66	!!----------------------------------------------------------------------
[2715]	67	!!
[2528]	68	INTEGER :: inum0 ! temprary units for 'mesh_mask.nc' file
	69	INTEGER :: inum1 ! temprary units for 'mesh.nc' file
	70	INTEGER :: inum2 ! temprary units for 'mask.nc' file
	71	INTEGER :: inum3 ! temprary units for 'mesh_hgr.nc' file
	72	INTEGER :: inum4 ! temprary units for 'mesh_zgr.nc' file
	73	CHARACTER(len=21) :: clnam0 ! filename (mesh and mask informations)
	74	CHARACTER(len=21) :: clnam1 ! filename (mesh informations)
	75	CHARACTER(len=21) :: clnam2 ! filename (mask informations)
	76	CHARACTER(len=21) :: clnam3 ! filename (horizontal mesh informations)
	77	CHARACTER(len=21) :: clnam4 ! filename (vertical mesh informations)
	78	INTEGER :: ji, jj, jk ! dummy loop indices
[3294]	79	! ! workspaces
	80	REAL(wp), POINTER, DIMENSION(:,: ) :: zprt, zprw
	81	REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepu, zdepv
[2715]	82	!!----------------------------------------------------------------------
[3294]	83	!
	84	IF( nn_timing == 1 ) CALL timing_start('dom_wri')
	85	!
	86	CALL wrk_alloc( jpi, jpj, zprt, zprw )
	87	CALL wrk_alloc( jpi, jpj, jpk, zdepu, zdepv )
	88	!
[1590]	89	IF(lwp) WRITE(numout,*)
	90	IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)'
	91	IF(lwp) WRITE(numout,*) '~~~~~~~'
	92
	93	clnam0 = 'mesh_mask' ! filename (mesh and mask informations)
	94	clnam1 = 'mesh' ! filename (mesh informations)
	95	clnam2 = 'mask' ! filename (mask informations)
	96	clnam3 = 'mesh_hgr' ! filename (horizontal mesh informations)
	97	clnam4 = 'mesh_zgr' ! filename (vertical mesh informations)
	98
	99	SELECT CASE ( MOD(nmsh, 3) )
[1161]	100	! ! ============================
	101	CASE ( 1 ) ! create 'mesh_mask.nc' file
	102	! ! ============================
	103	CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib )
	104	inum2 = inum0 ! put all the informations
	105	inum3 = inum0 ! in unit inum0
	106	inum4 = inum0
	107
	108	! ! ============================
	109	CASE ( 2 ) ! create 'mesh.nc' and
	110	! ! 'mask.nc' files
	111	! ! ============================
	112	CALL iom_open( TRIM(clnam1), inum1, ldwrt = .TRUE., kiolib = jprstlib )
	113	CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
	114	inum3 = inum1 ! put mesh informations
	115	inum4 = inum1 ! in unit inum1
	116	! ! ============================
[1590]	117	CASE ( 0 ) ! create 'mesh_hgr.nc'
[1161]	118	! ! 'mesh_zgr.nc' and
	119	! ! 'mask.nc' files
	120	! ! ============================
	121	CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
	122	CALL iom_open( TRIM(clnam3), inum3, ldwrt = .TRUE., kiolib = jprstlib )
	123	CALL iom_open( TRIM(clnam4), inum4, ldwrt = .TRUE., kiolib = jprstlib )
[2528]	124	!
[1161]	125	END SELECT
	126
	127	! ! masks (inum2)
	128	CALL iom_rstput( 0, 0, inum2, 'tmask', tmask, ktype = jp_i1 ) ! ! land-sea mask
	129	CALL iom_rstput( 0, 0, inum2, 'umask', umask, ktype = jp_i1 )
	130	CALL iom_rstput( 0, 0, inum2, 'vmask', vmask, ktype = jp_i1 )
	131	CALL iom_rstput( 0, 0, inum2, 'fmask', fmask, ktype = jp_i1 )
	132
[2715]	133	CALL dom_uniq( zprw, 'T' )
	134	zprt = tmask(:,:,1) * zprw ! ! unique point mask
[1161]	135	CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 )
[2715]	136	CALL dom_uniq( zprw, 'U' )
	137	zprt = umask(:,:,1) * zprw
[1161]	138	CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 )
[2715]	139	CALL dom_uniq( zprw, 'V' )
	140	zprt = vmask(:,:,1) * zprw
[1161]	141	CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 )
[2715]	142	CALL dom_uniq( zprw, 'F' )
	143	zprt = fmask(:,:,1) * zprw
[1161]	144	CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 )
	145
	146	! ! horizontal mesh (inum3)
	147	CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r4 ) ! ! latitude
	148	CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r4 )
	149	CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r4 )
	150	CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r4 )
	151
	152	CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r4 ) ! ! longitude
	153	CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r4 )
	154	CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r4 )
	155	CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r4 )
	156
	157	CALL iom_rstput( 0, 0, inum3, 'e1t', e1t, ktype = jp_r8 ) ! ! e1 scale factors
	158	CALL iom_rstput( 0, 0, inum3, 'e1u', e1u, ktype = jp_r8 )
	159	CALL iom_rstput( 0, 0, inum3, 'e1v', e1v, ktype = jp_r8 )
	160	CALL iom_rstput( 0, 0, inum3, 'e1f', e1f, ktype = jp_r8 )
	161
	162	CALL iom_rstput( 0, 0, inum3, 'e2t', e2t, ktype = jp_r8 ) ! ! e2 scale factors
	163	CALL iom_rstput( 0, 0, inum3, 'e2u', e2u, ktype = jp_r8 )
	164	CALL iom_rstput( 0, 0, inum3, 'e2v', e2v, ktype = jp_r8 )
	165	CALL iom_rstput( 0, 0, inum3, 'e2f', e2f, ktype = jp_r8 )
	166
	167	CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 ) ! ! coriolis factor
	168
[2528]	169	! note that mbkt is set to 1 over land ==> use surface tmask
	170	zprt(:,:) = tmask(:,:,1) * REAL( mbkt(:,:) , wp )
	171	CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 ) ! ! nb of ocean T-points
[1161]	172
	173	IF( ln_sco ) THEN ! s-coordinate
[3680]	174	CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )
	175	CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu )
[1161]	176	CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv )
	177	CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf )
[2528]	178	!
[1161]	179	CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt ) ! ! scaling coef.
	180	CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw )
	181	CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w )
	182	CALL iom_rstput( 0, 0, inum4, 'esigt', esigt )
	183	CALL iom_rstput( 0, 0, inum4, 'esigw', esigw )
[2528]	184	!
[4292]	185	CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 ) ! ! scale factors
	186	CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
	187	CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
	188	CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
[3680]	189	CALL iom_rstput( 0, 0, inum4, 'rx1', rx1 ) ! ! Max. grid stiffness ratio
[2528]	190	!
[4292]	191	CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d ) ! ! stretched system
	192	CALL iom_rstput( 0, 0, inum4, 'gdepw_1d' , gdepw_1d )
[1161]	193	ENDIF
	194
	195	IF( ln_zps ) THEN ! z-coordinate - partial steps
[2528]	196	!
[1590]	197	IF( nmsh <= 6 ) THEN ! ! 3D vertical scale factors
[4292]	198	CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )
	199	CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
	200	CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
	201	CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
[2528]	202	ELSE ! ! 2D masked bottom ocean scale factors
	203	DO jj = 1,jpj
	204	DO ji = 1,jpi
[4292]	205	e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
	206	e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
[2528]	207	END DO
	208	END DO
[1590]	209	CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp )
	210	CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp )
	211	END IF
[2528]	212	!
[1590]	213	IF( nmsh <= 3 ) THEN ! ! 3D depth
[4292]	214	CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r4 )
[2528]	215	DO jk = 1,jpk
	216	DO jj = 1, jpjm1
	217	DO ji = 1, fs_jpim1 ! vector opt.
[4292]	218	zdepu(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj ,jk) )
	219	zdepv(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji ,jj+1,jk) )
[2528]	220	END DO
	221	END DO
	222	END DO
[1590]	223	CALL lbc_lnk( zdepu, 'U', 1. ) ; CALL lbc_lnk( zdepv, 'V', 1. )
	224	CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 )
	225	CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 )
[4292]	226	CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r4 )
[1590]	227	ELSE ! ! 2D bottom depth
[2528]	228	DO jj = 1,jpj
	229	DO ji = 1,jpi
[4292]	230	zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj) ) * tmask(ji,jj,1)
	231	zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1)
[2528]	232	END DO
	233	END DO
	234	CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r4 )
	235	CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r4 )
[1590]	236	ENDIF
[2528]	237	!
[4292]	238	CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d ) ! ! reference z-coord.
	239	CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
	240	CALL iom_rstput( 0, 0, inum4, 'e3t_1d' , e3t_1d )
	241	CALL iom_rstput( 0, 0, inum4, 'e3w_1d' , e3w_1d )
[1161]	242	ENDIF
	243
	244	IF( ln_zco ) THEN
	245	! ! z-coordinate - full steps
[4292]	246	CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d ) ! ! depth
	247	CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
	248	CALL iom_rstput( 0, 0, inum4, 'e3t_1d' , e3t_1d ) ! ! scale factors
	249	CALL iom_rstput( 0, 0, inum4, 'e3w_1d' , e3w_1d )
[1161]	250	ENDIF
	251	! ! ============================
	252	! ! close the files
	253	! ! ============================
[1929]	254	SELECT CASE ( MOD(nmsh, 3) )
[1161]	255	CASE ( 1 )
	256	CALL iom_close( inum0 )
	257	CASE ( 2 )
	258	CALL iom_close( inum1 )
	259	CALL iom_close( inum2 )
[1929]	260	CASE ( 0 )
[1161]	261	CALL iom_close( inum2 )
	262	CALL iom_close( inum3 )
	263	CALL iom_close( inum4 )
	264	END SELECT
[2528]	265	!
[3294]	266	CALL wrk_dealloc( jpi, jpj, zprt, zprw )
	267	CALL wrk_dealloc( jpi, jpj, jpk, zdepu, zdepv )
[2715]	268	!
[3294]	269	IF( nn_timing == 1 ) CALL timing_stop('dom_wri')
	270	!
[1161]	271	END SUBROUTINE dom_wri
[3]	272
	273
[2715]	274	SUBROUTINE dom_uniq( puniq, cdgrd )
[1161]	275	!!----------------------------------------------------------------------
	276	!! * ROUTINE dom_uniq *
	277	!!
	278	!! ** Purpose : identify unique point of a grid (TUVF)
	279	!!
	280	!! ** Method : 1) aplly lbc_lnk on an array with different values for each element
	281	!! 2) check which elements have been changed
	282	!!----------------------------------------------------------------------
[2528]	283	!
[2715]	284	CHARACTER(len=1) , INTENT(in ) :: cdgrd !
	285	REAL(wp), DIMENSION(:,:), INTENT(inout) :: puniq !
	286	!
	287	REAL(wp) :: zshift ! shift value link to the process number
	288	INTEGER :: ji ! dummy loop indices
	289	LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) :: lldbl ! store whether each point is unique or not
[3294]	290	REAL(wp), POINTER, DIMENSION(:,:) :: ztstref
[1161]	291	!!----------------------------------------------------------------------
[3294]	292	!
	293	IF( nn_timing == 1 ) CALL timing_start('dom_uniq')
	294	!
	295	CALL wrk_alloc( jpi, jpj, ztstref )
	296	!
[1161]	297	! build an array with different values for each element
	298	! in mpp: make sure that these values are different even between process
	299	! -> apply a shift value according to the process number
	300	zshift = jpi * jpj * ( narea - 1 )
[2528]	301	ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
	302	!
[1161]	303	puniq(:,:) = ztstref(:,:) ! default definition
	304	CALL lbc_lnk( puniq, cdgrd, 1. ) ! apply boundary conditions
	305	lldbl(:,:,1) = puniq(:,:) == ztstref(:,:) ! check which values have been changed
[2528]	306	!
[1161]	307	puniq(:,:) = 1. ! default definition
	308	! fill only the inner part of the cpu with llbl converted into real
[2528]	309	puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
	310	!
[3294]	311	CALL wrk_dealloc( jpi, jpj, ztstref )
[2715]	312	!
[3294]	313	IF( nn_timing == 1 ) CALL timing_stop('dom_uniq')
	314	!
[2715]	315	END SUBROUTINE dom_uniq
[3]	316
	317	!!======================================================================
	318	END MODULE domwri

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

Download in other formats: