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.

domhgr.F90 in branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 @ 2380

Last change on this file since 2380 was 2380, checked in by acc, 13 years ago
nemo_v3_3beta. ORCA_R1 settings (step 2, see ticket #758). Introduces key_orca_r1 (46 level default, 75 level if key_orca_r1=75)
Property svn:keywords set to `Id`
File size: 31.2 KB

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

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

Download in other formats: