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 trunk/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90 @ 3648

Last change on this file since 3648 was 3294, checked in by rblod, 12 years ago
Merge of 3.4beta into the trunk
Property svn:keywords set to `Id`
File size: 16.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
[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
	174	CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt ) ! ! depth
	175	CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu )
	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	!
[1161]	185	CALL iom_rstput( 0, 0, inum4, 'e3t', e3t ) ! ! scale factors
	186	CALL iom_rstput( 0, 0, inum4, 'e3u', e3u )
	187	CALL iom_rstput( 0, 0, inum4, 'e3v', e3v )
	188	CALL iom_rstput( 0, 0, inum4, 'e3w', e3w )
[2528]	189	!
[1161]	190	CALL iom_rstput( 0, 0, inum4, 'gdept_0' , gdept_0 ) ! ! stretched system
	191	CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 )
	192	ENDIF
	193
	194	IF( ln_zps ) THEN ! z-coordinate - partial steps
[2528]	195	!
[1590]	196	IF( nmsh <= 6 ) THEN ! ! 3D vertical scale factors
	197	CALL iom_rstput( 0, 0, inum4, 'e3t', e3t )
	198	CALL iom_rstput( 0, 0, inum4, 'e3u', e3u )
	199	CALL iom_rstput( 0, 0, inum4, 'e3v', e3v )
	200	CALL iom_rstput( 0, 0, inum4, 'e3w', e3w )
[2528]	201	ELSE ! ! 2D masked bottom ocean scale factors
	202	DO jj = 1,jpj
	203	DO ji = 1,jpi
	204	e3tp(ji,jj) = e3t(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
	205	e3wp(ji,jj) = e3w(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)
	206	END DO
	207	END DO
[1590]	208	CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp )
	209	CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp )
	210	END IF
[2528]	211	!
[1590]	212	IF( nmsh <= 3 ) THEN ! ! 3D depth
	213	CALL iom_rstput( 0, 0, inum4, 'gdept', gdept, ktype = jp_r4 )
[2528]	214	DO jk = 1,jpk
	215	DO jj = 1, jpjm1
	216	DO ji = 1, fs_jpim1 ! vector opt.
	217	zdepu(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji+1,jj ,jk) )
	218	zdepv(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji ,jj+1,jk) )
	219	END DO
	220	END DO
	221	END DO
[1590]	222	CALL lbc_lnk( zdepu, 'U', 1. ) ; CALL lbc_lnk( zdepv, 'V', 1. )
	223	CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 )
	224	CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 )
	225	CALL iom_rstput( 0, 0, inum4, 'gdepw', gdepw, ktype = jp_r4 )
	226	ELSE ! ! 2D bottom depth
[2528]	227	DO jj = 1,jpj
	228	DO ji = 1,jpi
	229	zprt(ji,jj) = gdept(ji,jj,mbkt(ji,jj) ) * tmask(ji,jj,1)
	230	zprw(ji,jj) = gdepw(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1)
	231	END DO
	232	END DO
	233	CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r4 )
	234	CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r4 )
[1590]	235	ENDIF
[2528]	236	!
[1161]	237	CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0 ) ! ! reference z-coord.
	238	CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0 )
	239	CALL iom_rstput( 0, 0, inum4, 'e3t_0' , e3t_0 )
	240	CALL iom_rstput( 0, 0, inum4, 'e3w_0' , e3w_0 )
	241	ENDIF
	242
	243	IF( ln_zco ) THEN
	244	! ! z-coordinate - full steps
	245	CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0 ) ! ! depth
	246	CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0 )
	247	CALL iom_rstput( 0, 0, inum4, 'e3t_0' , e3t_0 ) ! ! scale factors
	248	CALL iom_rstput( 0, 0, inum4, 'e3w_0' , e3w_0 )
	249	ENDIF
	250	! ! ============================
	251	! ! close the files
	252	! ! ============================
[1929]	253	SELECT CASE ( MOD(nmsh, 3) )
[1161]	254	CASE ( 1 )
	255	CALL iom_close( inum0 )
	256	CASE ( 2 )
	257	CALL iom_close( inum1 )
	258	CALL iom_close( inum2 )
[1929]	259	CASE ( 0 )
[1161]	260	CALL iom_close( inum2 )
	261	CALL iom_close( inum3 )
	262	CALL iom_close( inum4 )
	263	END SELECT
[2528]	264	!
[3294]	265	CALL wrk_dealloc( jpi, jpj, zprt, zprw )
	266	CALL wrk_dealloc( jpi, jpj, jpk, zdepu, zdepv )
[2715]	267	!
[3294]	268	IF( nn_timing == 1 ) CALL timing_stop('dom_wri')
	269	!
[1161]	270	END SUBROUTINE dom_wri
[3]	271
	272
[2715]	273	SUBROUTINE dom_uniq( puniq, cdgrd )
[1161]	274	!!----------------------------------------------------------------------
	275	!! * ROUTINE dom_uniq *
	276	!!
	277	!! ** Purpose : identify unique point of a grid (TUVF)
	278	!!
	279	!! ** Method : 1) aplly lbc_lnk on an array with different values for each element
	280	!! 2) check which elements have been changed
	281	!!----------------------------------------------------------------------
[2528]	282	!
[2715]	283	CHARACTER(len=1) , INTENT(in ) :: cdgrd !
	284	REAL(wp), DIMENSION(:,:), INTENT(inout) :: puniq !
	285	!
	286	REAL(wp) :: zshift ! shift value link to the process number
	287	INTEGER :: ji ! dummy loop indices
	288	LOGICAL, DIMENSION(SIZE(puniq,1),SIZE(puniq,2),1) :: lldbl ! store whether each point is unique or not
[3294]	289	REAL(wp), POINTER, DIMENSION(:,:) :: ztstref
[1161]	290	!!----------------------------------------------------------------------
[3294]	291	!
	292	IF( nn_timing == 1 ) CALL timing_start('dom_uniq')
	293	!
	294	CALL wrk_alloc( jpi, jpj, ztstref )
	295	!
[1161]	296	! build an array with different values for each element
	297	! in mpp: make sure that these values are different even between process
	298	! -> apply a shift value according to the process number
	299	zshift = jpi * jpj * ( narea - 1 )
[2528]	300	ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
	301	!
[1161]	302	puniq(:,:) = ztstref(:,:) ! default definition
	303	CALL lbc_lnk( puniq, cdgrd, 1. ) ! apply boundary conditions
	304	lldbl(:,:,1) = puniq(:,:) == ztstref(:,:) ! check which values have been changed
[2528]	305	!
[1161]	306	puniq(:,:) = 1. ! default definition
	307	! fill only the inner part of the cpu with llbl converted into real
[2528]	308	puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp )
	309	!
[3294]	310	CALL wrk_dealloc( jpi, jpj, ztstref )
[2715]	311	!
[3294]	312	IF( nn_timing == 1 ) CALL timing_stop('dom_uniq')
	313	!
[2715]	314	END SUBROUTINE dom_uniq
[3]	315
	316	!!======================================================================
	317	END MODULE domwri

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

Download in other formats: