Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

domhgr.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago
The Dr Hook changes from my perl code.
File size: 34.8 KB

Rev	Line
[3]	1	MODULE domhgr
	2	!!==============================================================================
[93]	3	!! * MODULE domhgr *
[3]	4	!! Ocean initialization : domain initialization
	5	!!==============================================================================
[2528]	6	!! History : OPA ! 1988-03 (G. Madec) Original code
	7	!! 7.0 ! 1996-01 (G. Madec) terrain following coordinates
	8	!! 8.0 ! 1997-02 (G. Madec) print mesh informations
	9	!! 8.1 ! 1999-11 (M. Imbard) NetCDF format with IO-IPSL
	10	!! 8.2 ! 2000-08 (D. Ludicone) Reduced section at Bab el Mandeb
	11	!! - ! 2001-09 (M. Levy) eel config: grid in km, beta-plane
	12	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module, namelist
	13	!! - ! 2004-01 (A.M. Treguier, J.M. Molines) Case 4 (Mercator mesh)
	14	!! use of parameters in par_CONFIG-Rxx.h90, not in namelist
	15	!! - ! 2004-05 (A. Koch-Larrouy) Add Gyre configuration
[2715]	16	!! 4.0 ! 2011-02 (G. Madec) add cell surface (e1e2t)
[473]	17	!!----------------------------------------------------------------------
[3]	18
	19	!!----------------------------------------------------------------------
[2528]	20	!! dom_hgr : initialize the horizontal mesh
	21	!! hgr_read : read "coordinate" NetCDF file
[3]	22	!!----------------------------------------------------------------------
[2528]	23	USE dom_oce ! ocean space and time domain
	24	USE phycst ! physical constants
	25	USE in_out_manager ! I/O manager
	26	USE lib_mpp ! MPP library
[3294]	27	USE timing ! Timing
[3]	28
[11738]	29	USE yomhook, ONLY: lhook, dr_hook
	30	USE parkind1, ONLY: jprb, jpim
	31
[3]	32	IMPLICIT NONE
	33	PRIVATE
	34
[2528]	35	REAL(wp) :: glam0, gphi0 ! variables corresponding to parameters ppglam0 ppgphi0 set in par_oce
[3]	36
[2528]	37	PUBLIC dom_hgr ! called by domain.F90
	38
[3]	39	!!----------------------------------------------------------------------
[2715]	40	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
[6486]	41	!! $Id$
[2528]	42	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	43	!!----------------------------------------------------------------------
	44	CONTAINS
	45
	46	SUBROUTINE dom_hgr
	47	!!----------------------------------------------------------------------
	48	!! * ROUTINE dom_hgr *
	49	!!
	50	!! ** Purpose : Compute the geographical position (in degre) of the
	51	!! model grid-points, the horizontal scale factors (in meters) and
	52	!! the Coriolis factor (in s-1).
	53	!!
	54	!! ** Method : The geographical position of the model grid-points is
	55	!! defined from analytical functions, fslam and fsphi, the deriva-
	56	!! tives of which gives the horizontal scale factors e1,e2.
	57	!! Defining two function fslam and fsphi and their derivatives in
	58	!! the two horizontal directions (fse1 and fse2), the model grid-
	59	!! point position and scale factors are given by:
[81]	60	!! t-point:
	61	!! glamt(i,j) = fslam(i ,j ) e1t(i,j) = fse1(i ,j )
	62	!! gphit(i,j) = fsphi(i ,j ) e2t(i,j) = fse2(i ,j )
	63	!! u-point:
	64	!! glamu(i,j) = fslam(i+1/2,j ) e1u(i,j) = fse1(i+1/2,j )
	65	!! gphiu(i,j) = fsphi(i+1/2,j ) e2u(i,j) = fse2(i+1/2,j )
	66	!! v-point:
	67	!! glamv(i,j) = fslam(i ,j+1/2) e1v(i,j) = fse1(i ,j+1/2)
	68	!! gphiv(i,j) = fsphi(i ,j+1/2) e2v(i,j) = fse2(i ,j+1/2)
	69	!! f-point:
	70	!! glamf(i,j) = fslam(i+1/2,j+1/2) e1f(i,j) = fse1(i+1/2,j+1/2)
	71	!! gphif(i,j) = fsphi(i+1/2,j+1/2) e2f(i,j) = fse2(i+1/2,j+1/2)
[3]	72	!! Where fse1 and fse2 are defined by:
	73	!! fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2
	74	!! + di(fsphi) **2 )(i,j)
	75	!! fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2
	76	!! + dj(fsphi) **2 )(i,j)
	77	!!
	78	!! The coriolis factor is given at z-point by:
	79	!! ff = 2.omegasin(gphif) (in s-1)
	80	!!
	81	!! This routine is given as an example, it must be modified
	82	!! following the user s desiderata. nevertheless, the output as
	83	!! well as the way to compute the model grid-point position and
	84	!! horizontal scale factors must be respected in order to insure
	85	!! second order accuracy schemes.
	86	!!
	87	!! N.B. If the domain is periodic, verify that scale factors are also
	88	!! periodic, and the coriolis term again.
	89	!!
	90	!! ** Action : - define glamt, glamu, glamv, glamf: longitude of t-,
	91	!! u-, v- and f-points (in degre)
	92	!! - define gphit, gphiu, gphiv, gphit: latitude of t-,
	93	!! u-, v- and f-points (in degre)
	94	!! define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal
	95	!! scale factors (in meters) at t-, u-, v-, and f-points.
	96	!! define ff: coriolis factor at f-point
	97	!!
[473]	98	!! References : Marti, Madec and Delecluse, 1992, JGR
	99	!! Madec, Imbard, 1996, Clim. Dyn.
[3]	100	!!----------------------------------------------------------------------
[2528]	101	INTEGER :: ji, jj ! dummy loop indices
	102	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
	103	INTEGER :: ijeq ! index of equator T point (used in case 4)
	104	REAL(wp) :: zti, zui, zvi, zfi ! local scalars
	105	REAL(wp) :: ztj, zuj, zvj, zfj ! - -
	106	REAL(wp) :: zphi0, zbeta, znorme !
	107	REAL(wp) :: zarg, zf0, zminff, zmaxff
	108	REAL(wp) :: zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg
	109	REAL(wp) :: zphi1, zsin_alpha, zim05, zjm05
[5385]	110	INTEGER :: isrow ! index for ORCA1 starting row
[11738]	111	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
	112	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
	113	REAL(KIND=jprb) :: zhook_handle
[5385]	114
[11738]	115	CHARACTER(LEN=*), PARAMETER :: RoutineName='DOM_HGR'
	116
	117	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
	118
	119
[3]	120	!!----------------------------------------------------------------------
[3294]	121	!
	122	IF( nn_timing == 1 ) CALL timing_start('dom_hgr')
	123	!
[3]	124	IF(lwp) THEN
	125	WRITE(numout,*)
	126	WRITE(numout,*) 'dom_hgr : define the horizontal mesh from ithe following par_oce parameters '
	127	WRITE(numout,*) '~~~~~~~ type of horizontal mesh jphgr_msh = ', jphgr_msh
	128	WRITE(numout,*) ' position of the first row and ppglam0 = ', ppglam0
	129	WRITE(numout,*) ' column grid-point (degrees) ppgphi0 = ', ppgphi0
	130	WRITE(numout,*) ' zonal grid-spacing (degrees) ppe1_deg = ', ppe1_deg
	131	WRITE(numout,*) ' meridional grid-spacing (degrees) ppe2_deg = ', ppe2_deg
	132	WRITE(numout,*) ' zonal grid-spacing (meters) ppe1_m = ', ppe1_m
	133	WRITE(numout,*) ' meridional grid-spacing (meters) ppe2_m = ', ppe2_m
	134	ENDIF
	135
	136
	137	SELECT CASE( jphgr_msh ) ! type of horizontal mesh
	138
	139	CASE ( 0 ) ! curvilinear coordinate on the sphere read in coordinate.nc file
	140
	141	IF(lwp) WRITE(numout,*)
[81]	142	IF(lwp) WRITE(numout,*) ' curvilinear coordinate on the sphere read in "coordinate" file'
[239]	143
[81]	144	CALL hgr_read ! Defaultl option : NetCDF file
[3]	145
[81]	146	! ! =====================
	147	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
	148	! ! =====================
[2528]	149	IF( nn_cla == 0 ) THEN
[1273]	150	!
[236]	151	ii0 = 139 ; ii1 = 140 ! Gibraltar Strait (e2u = 20 km)
	152	ij0 = 102 ; ij1 = 102 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3
	153	IF(lwp) WRITE(numout,*)
[1273]	154	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar : e2u reduced to 20 km'
	155	!
[236]	156	ii0 = 160 ; ii1 = 160 ! Bab el Mandeb (e2u = 18 km)
[1273]	157	ij0 = 88 ; ij1 = 88 ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 18.e3
	158	e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3
[81]	159	IF(lwp) WRITE(numout,*)
[1273]	160	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb: e2u reduced to 30 km'
	161	IF(lwp) WRITE(numout,*) ' e1v reduced to 18 km'
[81]	162	ENDIF
[3]	163
[1707]	164	ii0 = 145 ; ii1 = 146 ! Danish Straits (e2u = 10 km)
	165	ij0 = 116 ; ij1 = 116 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e3
[81]	166	IF(lwp) WRITE(numout,*)
[1707]	167	IF(lwp) WRITE(numout,*) ' orca_r2: Danish Straits : e2u reduced to 10 km'
[81]	168	!
	169	ENDIF
	170
[2528]	171	! ! =====================
	172	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
	173	! ! =====================
[5385]	174	! This dirty section will be suppressed by simplification process: all this will come back in input files
[5506]	175	! Currently these hard-wired indices relate to configuration with
	176	! extend grid (jpjglo=332)
[5385]	177	! which had a grid-size of 362x292.
	178	!
[5506]	179	isrow = 332 - jpjglo
[5385]	180	!
	181	ii0 = 282 ; ii1 = 283 ! Gibraltar Strait (e2u = 20 km)
[6487]	182	ij0 = 241 - isrow ; ij1 = 241 - isrow ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3
[2528]	183	IF(lwp) WRITE(numout,*)
	184	IF(lwp) WRITE(numout,*) ' orca_r1: Gibraltar : e2u reduced to 20 km'
	185
[5385]	186	ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u = 10 km)
[6487]	187	ij0 = 248 - isrow ; ij1 = 248 - isrow ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e3
[2528]	188	IF(lwp) WRITE(numout,*)
	189	IF(lwp) WRITE(numout,*) ' orca_r1: Bhosporus : e2u reduced to 10 km'
	190
[5385]	191	ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v = 13 km)
[6487]	192	ij0 = 164 - isrow ; ij1 = 165 - isrow ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3
[2528]	193	IF(lwp) WRITE(numout,*)
	194	IF(lwp) WRITE(numout,*) ' orca_r1: Lombok : e1v reduced to 10 km'
	195
[5385]	196	ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on]
[6487]	197	ij0 = 164 - isrow ; ij1 = 165 - isrow ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 8.e3
[2528]	198	IF(lwp) WRITE(numout,*)
	199	IF(lwp) WRITE(numout,*) ' orca_r1: Sumba : e1v reduced to 8 km'
	200
[5385]	201	ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v = 13 km)
[6487]	202	ij0 = 164 - isrow ; ij1 = 165 - isrow ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3
[2528]	203	IF(lwp) WRITE(numout,*)
	204	IF(lwp) WRITE(numout,*) ' orca_r1: Ombai : e1v reduced to 13 km'
	205
[5385]	206	ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v = 20 km)
[6487]	207	ij0 = 164 - isrow ; ij1 = 145 - isrow ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3
[2528]	208	IF(lwp) WRITE(numout,*)
	209	IF(lwp) WRITE(numout,*) ' orca_r1: Timor Passage : e1v reduced to 20 km'
	210
[5385]	211	ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v = 30 km)
[6487]	212	ij0 = 181 - isrow ; ij1 = 182 - isrow ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3
[2528]	213	IF(lwp) WRITE(numout,*)
	214	IF(lwp) WRITE(numout,*) ' orca_r1: W Halmahera : e1v reduced to 30 km'
	215
[5385]	216	ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v = 50 km)
[6487]	217	ij0 = 181 - isrow ; ij1 = 182 - isrow ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3
[2528]	218	IF(lwp) WRITE(numout,*)
	219	IF(lwp) WRITE(numout,*) ' orca_r1: E Halmahera : e1v reduced to 50 km'
	220	!
	221	!
	222	ENDIF
	223
[81]	224	! ! ======================
	225	IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN ! ORCA R05 configuration
	226	! ! ======================
	227	ii0 = 563 ; ii1 = 564 ! Gibraltar Strait (e2u = 20 km)
	228	ij0 = 327 ; ij1 = 327 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3
	229	IF(lwp) WRITE(numout,*)
	230	IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e2u at the Gibraltar Strait'
	231	!
[473]	232	ii0 = 627 ; ii1 = 628 ! Bosphore Strait (e2u = 10 km)
	233	ij0 = 343 ; ij1 = 343 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e3
	234	IF(lwp) WRITE(numout,*)
	235	IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e2u at the Bosphore Strait'
	236	!
	237	ii0 = 93 ; ii1 = 94 ! Sumba Strait (e2u = 40 km)
	238	ij0 = 232 ; ij1 = 232 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 40.e3
	239	IF(lwp) WRITE(numout,*)
	240	IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e2u at the Sumba Strait'
	241	!
	242	ii0 = 103 ; ii1 = 103 ! Ombai Strait (e2u = 15 km)
	243	ij0 = 232 ; ij1 = 232 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 15.e3
	244	IF(lwp) WRITE(numout,*)
	245	IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e2u at the Ombai Strait'
	246	!
	247	ii0 = 15 ; ii1 = 15 ! Palk Strait (e2u = 10 km)
	248	ij0 = 270 ; ij1 = 270 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e3
	249	IF(lwp) WRITE(numout,*)
	250	IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e2u at the Palk Strait'
	251	!
	252	ii0 = 87 ; ii1 = 87 ! Lombok Strait (e1v = 10 km)
	253	ij0 = 232 ; ij1 = 233 ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e3
	254	IF(lwp) WRITE(numout,*)
	255	IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e1v at the Lombok Strait'
	256	!
	257	!
	258	ii0 = 662 ; ii1 = 662 ! Bab el Mandeb (e1v = 25 km)
	259	ij0 = 276 ; ij1 = 276 ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 25.e3
	260	IF(lwp) WRITE(numout,*)
	261	IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e1v at the Bab el Mandeb'
	262	!
[81]	263	ENDIF
	264
	265
[3]	266	! N.B. : General case, lat and long function of both i and j indices:
	267	! e1t(ji,jj) = ra * rad * SQRT( ( cos( radgphit(ji,jj) ) fsdila( zti, ztj ) )**2 &
	268	! + ( fsdiph( zti, ztj ) )**2 )
	269	! e1u(ji,jj) = ra * rad * SQRT( ( cos( radgphiu(ji,jj) ) fsdila( zui, zuj ) )**2 &
	270	! + ( fsdiph( zui, zuj ) )**2 )
	271	! e1v(ji,jj) = ra * rad * SQRT( ( cos( radgphiv(ji,jj) ) fsdila( zvi, zvj ) )**2 &
	272	! + ( fsdiph( zvi, zvj ) )**2 )
	273	! e1f(ji,jj) = ra * rad * SQRT( ( cos( radgphif(ji,jj) ) fsdila( zfi, zfj ) )**2 &
	274	! + ( fsdiph( zfi, zfj ) )**2 )
	275	!
	276	! e2t(ji,jj) = ra * rad * SQRT( ( cos( radgphit(ji,jj) ) fsdjla( zti, ztj ) )**2 &
	277	! + ( fsdjph( zti, ztj ) )**2 )
	278	! e2u(ji,jj) = ra * rad * SQRT( ( cos( radgphiu(ji,jj) ) fsdjla( zui, zuj ) )**2 &
	279	! + ( fsdjph( zui, zuj ) )**2 )
	280	! e2v(ji,jj) = ra * rad * SQRT( ( cos( radgphiv(ji,jj) ) fsdjla( zvi, zvj ) )**2 &
	281	! + ( fsdjph( zvi, zvj ) )**2 )
	282	! e2f(ji,jj) = ra * rad * SQRT( ( cos( radgphif(ji,jj) ) fsdjla( zfi, zfj ) )**2 &
	283	! + ( fsdjph( zfi, zfj ) )**2 )
	284
	285
	286	CASE ( 1 ) ! geographical mesh on the sphere with regular grid-spacing
	287
	288	IF(lwp) WRITE(numout,*)
	289	IF(lwp) WRITE(numout,*) ' geographical mesh on the sphere with regular grid-spacing'
	290	IF(lwp) WRITE(numout,*) ' given by ppe1_deg and ppe2_deg'
	291
	292	DO jj = 1, jpj
	293	DO ji = 1, jpi
	294	zti = FLOAT( ji - 1 + nimpp - 1 ) ; ztj = FLOAT( jj - 1 + njmpp - 1 )
	295	zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5 ; zuj = FLOAT( jj - 1 + njmpp - 1 )
	296	zvi = FLOAT( ji - 1 + nimpp - 1 ) ; zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
	297	zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5 ; zfj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5
	298	! Longitude
	299	glamt(ji,jj) = ppglam0 + ppe1_deg * zti
	300	glamu(ji,jj) = ppglam0 + ppe1_deg * zui
	301	glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
	302	glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
	303	! Latitude
	304	gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj
	305	gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj
	306	gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj
	307	gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj
	308	! e1
	309	e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
	310	e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
	311	e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
	312	e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
	313	! e2
	314	e2t(ji,jj) = ra * rad * ppe2_deg
	315	e2u(ji,jj) = ra * rad * ppe2_deg
	316	e2v(ji,jj) = ra * rad * ppe2_deg
	317	e2f(ji,jj) = ra * rad * ppe2_deg
	318	END DO
	319	END DO
	320
	321
	322	CASE ( 2:3 ) ! f- or beta-plane with regular grid-spacing
	323
	324	IF(lwp) WRITE(numout,*)
	325	IF(lwp) WRITE(numout,*) ' f- or beta-plane with regular grid-spacing'
	326	IF(lwp) WRITE(numout,*) ' given by ppe1_m and ppe2_m'
	327
	328	! Position coordinates (in kilometers)
	329	! ==========
	330	glam0 = 0.e0
	331	gphi0 = - ppe2_m * 1.e-3
[389]	332
[4147]	333	#if defined key_agrif
	334	IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN ! for EEL6 configuration only
	335	IF( .NOT. Agrif_Root() ) THEN
	336	glam0 = Agrif_Parent(glam0) + (Agrif_ix())Agrif_Parent(ppe1_m) 1.e-3
	337	gphi0 = Agrif_Parent(gphi0) + (Agrif_iy())Agrif_Parent(ppe2_m) 1.e-3
	338	ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox()
	339	ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()
	340	ENDIF
[389]	341	ENDIF
	342	#endif
[3]	343	DO jj = 1, jpj
	344	DO ji = 1, jpi
	345	glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) )
	346	glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 )
	347	glamv(ji,jj) = glamt(ji,jj)
	348	glamf(ji,jj) = glamu(ji,jj)
	349
	350	gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) )
	351	gphiu(ji,jj) = gphit(ji,jj)
	352	gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 )
	353	gphif(ji,jj) = gphiv(ji,jj)
	354	END DO
	355	END DO
	356
	357	! Horizontal scale factors (in meters)
	358	! ======
	359	e1t(:,:) = ppe1_m ; e2t(:,:) = ppe2_m
	360	e1u(:,:) = ppe1_m ; e2u(:,:) = ppe2_m
	361	e1v(:,:) = ppe1_m ; e2v(:,:) = ppe2_m
	362	e1f(:,:) = ppe1_m ; e2f(:,:) = ppe2_m
	363
	364	CASE ( 4 ) ! geographical mesh on the sphere, isotropic MERCATOR type
	365
	366	IF(lwp) WRITE(numout,*)
	367	IF(lwp) WRITE(numout,*) ' geographical mesh on the sphere, MERCATOR type'
	368	IF(lwp) WRITE(numout,*) ' longitudinal/latitudinal spacing given by ppe1_deg'
[473]	369	IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
[3]	370
	371	! Find index corresponding to the equator, given the grid spacing e1_deg
	372	! and the (approximate) southern latitude ppgphi0.
	373	! This way we ensure that the equator is at a "T / U" point, when in the domain.
	374	! The formula should work even if the equator is outside the domain.
	375	zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
[29]	376	ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
[224]	377	IF( ppgphi0 > 0 ) ijeq = -ijeq
[3]	378
[29]	379	IF(lwp) WRITE(numout,*) ' Index of the equator on the MERCATOR grid:', ijeq
[3]	380
	381	DO jj = 1, jpj
	382	DO ji = 1, jpi
[29]	383	zti = FLOAT( ji - 1 + nimpp - 1 ) ; ztj = FLOAT( jj - ijeq + njmpp - 1 )
	384	zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5 ; zuj = FLOAT( jj - ijeq + njmpp - 1 )
	385	zvi = FLOAT( ji - 1 + nimpp - 1 ) ; zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
	386	zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5 ; zfj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5
[3]	387	! Longitude
	388	glamt(ji,jj) = ppglam0 + ppe1_deg * zti
	389	glamu(ji,jj) = ppglam0 + ppe1_deg * zui
	390	glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
	391	glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
	392	! Latitude
	393	gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg rad ztj ) )
[93]	394	gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg rad zuj ) )
	395	gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg rad zvj ) )
	396	gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg rad zfj ) )
[3]	397	! e1
	398	e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
	399	e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
	400	e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
	401	e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
	402	! e2
	403	e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
	404	e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
	405	e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
	406	e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
	407	END DO
	408	END DO
	409
[93]	410	CASE ( 5 ) ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
	411
	412	IF(lwp) WRITE(numout,*)
	413	IF(lwp) WRITE(numout,*) ' beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
	414	IF(lwp) WRITE(numout,*) ' given by ppe1_m and ppe2_m'
	415
	416	! Position coordinates (in kilometers)
	417	! ==========
	418
[167]	419	! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
[93]	420	zlam1 = -85
	421	zphi1 = 29
[167]	422	! resolution in meters
	423	ze1 = 106000. / FLOAT(jp_cfg)
	424	! benchmark: forced the resolution to be about 100 km
	425	IF( nbench /= 0 ) ze1 = 106000.e0
[93]	426	zsin_alpha = - SQRT( 2. ) / 2.
	427	zcos_alpha = SQRT( 2. ) / 2.
	428	ze1deg = ze1 / (ra * rad)
	429	IF( nbench /= 0 ) ze1deg = ze1deg / FLOAT(jp_cfg) ! benchmark: keep the lat/+lon
[167]	430	! ! at the right jp_cfg resolution
[93]	431	glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjglo-2 )
	432	gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjglo-2 )
	433
[516]	434	IF( nprint==1 .AND. lwp ) THEN
	435	WRITE(numout,*) ' ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
	436	WRITE(numout,*) ' ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
	437	ENDIF
[93]	438
	439	DO jj = 1, jpj
	440	DO ji = 1, jpi
	441	zim1 = FLOAT( ji + nimpp - 1 ) - 1. ; zim05 = FLOAT( ji + nimpp - 1 ) - 1.5
	442	zjm1 = FLOAT( jj + njmpp - 1 ) - 1. ; zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5
	443
	444	glamf(ji,jj) = glam0 + zim1 * ze1deg * zcos_alpha + zjm1 * ze1deg * zsin_alpha
	445	gphif(ji,jj) = gphi0 - zim1 * ze1deg * zsin_alpha + zjm1 * ze1deg * zcos_alpha
	446
	447	glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
	448	gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
	449
	450	glamu(ji,jj) = glam0 + zim1 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
	451	gphiu(ji,jj) = gphi0 - zim1 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
	452
	453	glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1 * ze1deg * zsin_alpha
	454	gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1 * ze1deg * zcos_alpha
	455	END DO
	456	END DO
	457
	458	! Horizontal scale factors (in meters)
	459	! ======
	460	e1t(:,:) = ze1 ; e2t(:,:) = ze1
	461	e1u(:,:) = ze1 ; e2u(:,:) = ze1
	462	e1v(:,:) = ze1 ; e2v(:,:) = ze1
	463	e1f(:,:) = ze1 ; e2f(:,:) = ze1
	464
[3]	465	CASE DEFAULT
[473]	466	WRITE(ctmp1,*) ' bad flag value for jphgr_msh = ', jphgr_msh
	467	CALL ctl_stop( ctmp1 )
[3]	468
	469	END SELECT
[2715]	470
	471	! T-cell surface
	472	! --------------
	473	e1e2t(:,:) = e1t(:,:) * e2t(:,:)
[4366]	474
	475	! Useful shortcuts (JC: note the duplicated e2e2t array ! Need some cleaning)
	476	! ---------------------------------------------------------------------------
	477	e12t (:,:) = e1t(:,:) * e2t(:,:)
	478	e12u (:,:) = e1u(:,:) * e2u(:,:)
	479	e12v (:,:) = e1v(:,:) * e2v(:,:)
	480	e12f (:,:) = e1f(:,:) * e2f(:,:)
	481	r1_e12t (:,:) = 1._wp / e12t(:,:)
	482	r1_e12u (:,:) = 1._wp / e12u(:,:)
	483	r1_e12v (:,:) = 1._wp / e12v(:,:)
	484	r1_e12f (:,:) = 1._wp / e12f(:,:)
	485	re2u_e1u(:,:) = e2u(:,:) / e1u(:,:)
	486	re1v_e2v(:,:) = e1v(:,:) / e2v(:,:)
[5123]	487	r1_e1t (:,:) = 1._wp / e1t(:,:)
	488	r1_e1u (:,:) = 1._wp / e1u(:,:)
	489	r1_e1v (:,:) = 1._wp / e1v(:,:)
	490	r1_e1f (:,:) = 1._wp / e1f(:,:)
	491	r1_e2t (:,:) = 1._wp / e2t(:,:)
	492	r1_e2u (:,:) = 1._wp / e2u(:,:)
	493	r1_e2v (:,:) = 1._wp / e2v(:,:)
	494	r1_e2f (:,:) = 1._wp / e2f(:,:)
[3]	495
	496	! Control printing : Grid informations (if not restart)
	497	! ----------------
	498
[516]	499	IF( lwp .AND. .NOT.ln_rstart ) THEN
[3]	500	WRITE(numout,*)
	501	WRITE(numout,*) ' longitude and e1 scale factors'
	502	WRITE(numout,*) ' ------------------------------'
	503	WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1), &
	504	glamv(ji,1), glamf(ji,1), &
	505	e1t(ji,1), e1u(ji,1), &
	506	e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
	507	9300 FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x, &
	508	f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
	509
	510	WRITE(numout,*)
	511	WRITE(numout,*) ' latitude and e2 scale factors'
	512	WRITE(numout,*) ' -----------------------------'
	513	WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj), &
	514	& gphiv(1,jj), gphif(1,jj), &
	515	& e2t (1,jj), e2u (1,jj), &
	516	& e2v (1,jj), e2f (1,jj), jj = 1, jpj, 10 )
	517	ENDIF
	518
	519
	520	IF( nprint == 1 .AND. lwp ) THEN
	521	WRITE(numout,*) ' e1u e2u '
	522	CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
	523	CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
	524	WRITE(numout,*) ' e1v e2v '
	525	CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
	526	CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
	527	WRITE(numout,*) ' e1f e2f '
	528	CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
	529	CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
	530	ENDIF
	531
	532
	533	! ================= !
	534	! Coriolis factor !
	535	! ================= !
	536
	537	SELECT CASE( jphgr_msh ) ! type of horizontal mesh
	538
	539	CASE ( 0, 1, 4 ) ! mesh on the sphere
	540
	541	ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) )
	542
	543	CASE ( 2 ) ! f-plane at ppgphi0
	544
	545	ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
	546
	547	IF(lwp) WRITE(numout,*) ' f-plane: Coriolis parameter = constant = ', ff(1,1)
	548
	549	CASE ( 3 ) ! beta-plane
	550
[187]	551	zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra ! beta at latitude ppgphi0
	552	zphi0 = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad ) ! latitude of the first row F-points
[389]	553
[4147]	554	#if defined key_agrif
	555	IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN ! for EEL6 configuration only
	556	IF( .NOT. Agrif_Root() ) THEN
[6487]	557	zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) &
	558	& / (ra * rad)
[4147]	559	ENDIF
[389]	560	ENDIF
	561	#endif
[187]	562	zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south
[3]	563
[187]	564	ff(:,:) = ( zf0 + zbeta * gphif(:,:) * 1.e+3 ) ! f = f0 +beta* y ( y=0 at south)
[389]	565
[516]	566	IF(lwp) THEN
	567	WRITE(numout,*)
	568	WRITE(numout,*) ' Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
	569	WRITE(numout,*) ' Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
	570	ENDIF
[434]	571	IF( lk_mpp ) THEN
	572	zminff=ff(nldi,nldj)
	573	zmaxff=ff(nldi,nlej)
	574	CALL mpp_min( zminff ) ! min over the global domain
	575	CALL mpp_max( zmaxff ) ! max over the global domain
[516]	576	IF(lwp) WRITE(numout,*) ' Coriolis parameter varies globally from ', zminff,' to ', zmaxff
[434]	577	END IF
[3]	578
[434]	579	CASE ( 5 ) ! beta-plane and rotated domain (gyre configuration)
[93]	580
	581	zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra ! beta at latitude ppgphi0
	582	zphi0 = 15.e0 ! latitude of the first row F-points
	583	zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south
	584
	585	ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south)
	586
[516]	587	IF(lwp) THEN
	588	WRITE(numout,*)
	589	WRITE(numout,*) ' Beta-plane and rotated domain : '
	590	WRITE(numout,*) ' Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
	591	ENDIF
	592
[434]	593	IF( lk_mpp ) THEN
	594	zminff=ff(nldi,nldj)
	595	zmaxff=ff(nldi,nlej)
	596	CALL mpp_min( zminff ) ! min over the global domain
	597	CALL mpp_max( zmaxff ) ! max over the global domain
[516]	598	IF(lwp) WRITE(numout,*) ' Coriolis parameter varies globally from ', zminff,' to ', zmaxff
[434]	599	END IF
[93]	600
[3]	601	END SELECT
	602
	603
	604	! Control of domain for symetrical condition
	605	! ------------------------------------------
	606	! The equator line must be the latitude coordinate axe
	607
	608	IF( nperio == 2 ) THEN
	609	znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )
[473]	610	IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
[3]	611	ENDIF
[3294]	612	!
	613	IF( nn_timing == 1 ) CALL timing_stop('dom_hgr')
	614	!
[11738]	615	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
[3]	616	END SUBROUTINE dom_hgr
	617
	618
[81]	619	SUBROUTINE hgr_read
[3]	620	!!---------------------------------------------------------------------
[81]	621	!! * ROUTINE hgr_read *
[3]	622	!!
	623	!! ** Purpose : Read a coordinate file in NetCDF format
	624	!!
	625	!! ** Method : The mesh file has been defined trough a analytical
	626	!! or semi-analytical method. It is read in a NetCDF file.
	627	!!
	628	!!----------------------------------------------------------------------
[473]	629	USE iom
[3]	630
[473]	631	INTEGER :: inum ! temporary logical unit
[11738]	632	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
	633	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
	634	REAL(KIND=jprb) :: zhook_handle
	635
	636	CHARACTER(LEN=*), PARAMETER :: RoutineName='HGR_READ'
	637
	638	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
	639
[3]	640	!!----------------------------------------------------------------------
	641
	642	IF(lwp) THEN
	643	WRITE(numout,*)
[81]	644	WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
[473]	645	WRITE(numout,*) '~~~~~~~~ jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
[3]	646	ENDIF
[473]	647
	648	CALL iom_open( 'coordinates', inum )
	649
[5118]	650	CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr )
	651	CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr )
	652	CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr )
	653	CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr )
[473]	654
[5118]	655	CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr )
	656	CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr )
	657	CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr )
	658	CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr )
[473]	659
[5118]	660	CALL iom_get( inum, jpdom_data, 'e1t', e1t, lrowattr=ln_use_jattr )
	661	CALL iom_get( inum, jpdom_data, 'e1u', e1u, lrowattr=ln_use_jattr )
	662	CALL iom_get( inum, jpdom_data, 'e1v', e1v, lrowattr=ln_use_jattr )
	663	CALL iom_get( inum, jpdom_data, 'e1f', e1f, lrowattr=ln_use_jattr )
[473]	664
[5118]	665	CALL iom_get( inum, jpdom_data, 'e2t', e2t, lrowattr=ln_use_jattr )
	666	CALL iom_get( inum, jpdom_data, 'e2u', e2u, lrowattr=ln_use_jattr )
	667	CALL iom_get( inum, jpdom_data, 'e2v', e2v, lrowattr=ln_use_jattr )
	668	CALL iom_get( inum, jpdom_data, 'e2f', e2f, lrowattr=ln_use_jattr )
[473]	669
	670	CALL iom_close( inum )
	671
[4990]	672	! need to be define for the extended grid south of -80S
	673	! some point are undefined but you need to have e1 and e2 .NE. 0
	674	WHERE (e1t==0.0_wp)
	675	e1t=1.0e2
	676	END WHERE
	677	WHERE (e1v==0.0_wp)
	678	e1v=1.0e2
	679	END WHERE
	680	WHERE (e1u==0.0_wp)
	681	e1u=1.0e2
	682	END WHERE
	683	WHERE (e1f==0.0_wp)
	684	e1f=1.0e2
	685	END WHERE
	686	WHERE (e2t==0.0_wp)
	687	e2t=1.0e2
	688	END WHERE
	689	WHERE (e2v==0.0_wp)
	690	e2v=1.0e2
	691	END WHERE
	692	WHERE (e2u==0.0_wp)
	693	e2u=1.0e2
	694	END WHERE
	695	WHERE (e2f==0.0_wp)
	696	e2f=1.0e2
	697	END WHERE
	698
[11738]	699	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
[473]	700	END SUBROUTINE hgr_read
	701
[3]	702	!!======================================================================
	703	END MODULE domhgr

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 11738

Download in other formats: