Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

geo2ocean.F90

Last change on this file was 14433, checked in by smasson, 3 years ago
trunk: merge dev_r14312_MPI_Interface into the trunk, #2598
Property svn:keywords set to `Id`
File size: 21.9 KB

Rev	Line
[3]	1	MODULE geo2ocean
	2	!!======================================================================
	3	!! * MODULE geo2ocean *
[144]	4	!! Ocean mesh : ???
[1218]	5	!!======================================================================
	6	!! History : OPA ! 07-1996 (O. Marti) Original code
[6140]	7	!! NEMO 1.0 ! 06-2006 (G. Madec ) Free form, F90 + opt.
	8	!! ! 04-2007 (S. Masson) angle: Add T, F points and bugfix in cos lateral boundary
	9	!! 3.0 ! 07-2008 (G. Madec) geo2oce suppress lon/lat agruments
	10	!! 3.7 ! 11-2015 (G. Madec) remove the unused repere and repcmo routines
[1218]	11	!!----------------------------------------------------------------------
[3]	12	!!----------------------------------------------------------------------
[6140]	13	!! rot_rep : Rotate the Repere: geographic grid <==> stretched coordinates grid
	14	!! angle :
	15	!! geo2oce :
	16	!! oce2geo :
[3]	17	!!----------------------------------------------------------------------
[6140]	18	USE dom_oce ! mesh and scale factors
	19	USE phycst ! physical constants
	20	!
	21	USE in_out_manager ! I/O manager
	22	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	23	USE lib_mpp ! MPP library
[3]	24
	25	IMPLICIT NONE
[1218]	26	PRIVATE
[3]	27
[6140]	28	PUBLIC rot_rep ! called in sbccpl, fldread, and cyclone
	29	PUBLIC geo2oce ! called in sbccpl
	30	PUBLIC oce2geo ! called in sbccpl
	31	PUBLIC obs_rot ! called in obs_rot_vel and obs_write
[3]	32
[6140]	33	! ! cos/sin between model grid lines and NP direction
	34	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: gsint, gcost ! at T point
	35	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: gsinu, gcosu ! at U point
	36	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: gsinv, gcosv ! at V point
	37	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: gsinf, gcosf ! at F point
[2528]	38
[2715]	39	LOGICAL , SAVE, DIMENSION(4) :: linit = .FALSE.
	40	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gsinlon, gcoslon, gsinlat, gcoslat
	41
[6140]	42	LOGICAL :: lmust_init = .TRUE. !: used to initialize the cos/sin variables (see above)
[672]	43
[1218]	44	!! * Substitutions
[12377]	45	# include "do_loop_substitute.h90"
[1218]	46	!!----------------------------------------------------------------------
[9598]	47	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1226]	48	!! $Id$
[10068]	49	!! Software governed by the CeCILL license (see ./LICENSE)
[1218]	50	!!----------------------------------------------------------------------
[3]	51	CONTAINS
	52
[685]	53	SUBROUTINE rot_rep ( pxin, pyin, cd_type, cdtodo, prot )
[672]	54	!!----------------------------------------------------------------------
	55	!! * ROUTINE rot_rep *
	56	!!
	57	!! ** Purpose : Rotate the Repere: Change vector componantes between
	58	!! geographic grid <--> stretched coordinates grid.
	59	!!----------------------------------------------------------------------
[6140]	60	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pxin, pyin ! vector componantes
	61	CHARACTER(len=1), INTENT(in ) :: cd_type ! define the nature of pt2d array grid-points
	62	CHARACTER(len=5), INTENT(in ) :: cdtodo ! type of transpormation:
	63	! ! 'en->i' = east-north to i-component
	64	! ! 'en->j' = east-north to j-component
	65	! ! 'ij->e' = (i,j) components to east
	66	! ! 'ij->n' = (i,j) components to north
	67	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: prot
[672]	68	!!----------------------------------------------------------------------
[6140]	69	!
	70	IF( lmust_init ) THEN ! at 1st call only: set gsin. & gcos.
[3]	71	IF(lwp) WRITE(numout,*)
[6140]	72	IF(lwp) WRITE(numout,*) ' rot_rep: coordinate transformation : geographic <==> model (i,j)-components'
	73	IF(lwp) WRITE(numout,*) ' ~~~~~~~~ '
[2715]	74	!
[10370]	75	CALL angle( glamt, gphit, glamu, gphiu, glamv, gphiv, glamf, gphif ) ! initialization of the transformation
[672]	76	lmust_init = .FALSE.
[3]	77	ENDIF
[6140]	78	!
	79	SELECT CASE( cdtodo ) ! type of rotation
	80	!
	81	CASE( 'en->i' ) ! east-north to i-component
[672]	82	SELECT CASE (cd_type)
[7753]	83	CASE ('T') ; prot(:,:) = pxin(:,:) * gcost(:,:) + pyin(:,:) * gsint(:,:)
	84	CASE ('U') ; prot(:,:) = pxin(:,:) * gcosu(:,:) + pyin(:,:) * gsinu(:,:)
	85	CASE ('V') ; prot(:,:) = pxin(:,:) * gcosv(:,:) + pyin(:,:) * gsinv(:,:)
	86	CASE ('F') ; prot(:,:) = pxin(:,:) * gcosf(:,:) + pyin(:,:) * gsinf(:,:)
[672]	87	CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
	88	END SELECT
[6140]	89	CASE ('en->j') ! east-north to j-component
[672]	90	SELECT CASE (cd_type)
[7753]	91	CASE ('T') ; prot(:,:) = pyin(:,:) * gcost(:,:) - pxin(:,:) * gsint(:,:)
	92	CASE ('U') ; prot(:,:) = pyin(:,:) * gcosu(:,:) - pxin(:,:) * gsinu(:,:)
	93	CASE ('V') ; prot(:,:) = pyin(:,:) * gcosv(:,:) - pxin(:,:) * gsinv(:,:)
	94	CASE ('F') ; prot(:,:) = pyin(:,:) * gcosf(:,:) - pxin(:,:) * gsinf(:,:)
[672]	95	CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
	96	END SELECT
[6140]	97	CASE ('ij->e') ! (i,j)-components to east
[672]	98	SELECT CASE (cd_type)
[7753]	99	CASE ('T') ; prot(:,:) = pxin(:,:) * gcost(:,:) - pyin(:,:) * gsint(:,:)
	100	CASE ('U') ; prot(:,:) = pxin(:,:) * gcosu(:,:) - pyin(:,:) * gsinu(:,:)
	101	CASE ('V') ; prot(:,:) = pxin(:,:) * gcosv(:,:) - pyin(:,:) * gsinv(:,:)
	102	CASE ('F') ; prot(:,:) = pxin(:,:) * gcosf(:,:) - pyin(:,:) * gsinf(:,:)
[672]	103	CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
	104	END SELECT
[6140]	105	CASE ('ij->n') ! (i,j)-components to north
[672]	106	SELECT CASE (cd_type)
[7753]	107	CASE ('T') ; prot(:,:) = pyin(:,:) * gcost(:,:) + pxin(:,:) * gsint(:,:)
	108	CASE ('U') ; prot(:,:) = pyin(:,:) * gcosu(:,:) + pxin(:,:) * gsinu(:,:)
	109	CASE ('V') ; prot(:,:) = pyin(:,:) * gcosv(:,:) + pxin(:,:) * gsinv(:,:)
	110	CASE ('F') ; prot(:,:) = pyin(:,:) * gcosf(:,:) + pxin(:,:) * gsinf(:,:)
[672]	111	CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
	112	END SELECT
	113	CASE DEFAULT ; CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' )
[6140]	114	!
[672]	115	END SELECT
[6140]	116	!
[685]	117	END SUBROUTINE rot_rep
[3]	118
	119
[10370]	120	SUBROUTINE angle( plamt, pphit, plamu, pphiu, plamv, pphiv, plamf, pphif )
[3]	121	!!----------------------------------------------------------------------
	122	!! * ROUTINE angle *
	123	!!
[672]	124	!! ** Purpose : Compute angles between model grid lines and the North direction
[3]	125	!!
[6140]	126	!! ** Method : sinus and cosinus of the angle between the north-south axe
	127	!! and the j-direction at t, u, v and f-points
	128	!! dot and cross products are used to obtain cos and sin, resp.
[3]	129	!!
[6140]	130	!! ** Action : - gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf
[3]	131	!!----------------------------------------------------------------------
[10370]	132	! WARNING: for an unexplained reason, we need to pass all glam, gphi arrays as input parameters in
	133	! order to get AGRIF working with -03 compilation option
	134	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: plamt, pphit, plamu, pphiu, plamv, pphiv, plamf, pphif
	135	!
[6140]	136	INTEGER :: ji, jj ! dummy loop indices
	137	INTEGER :: ierr ! local integer
	138	REAL(wp) :: zlam, zphi ! local scalars
	139	REAL(wp) :: zlan, zphh ! - -
	140	REAL(wp) :: zxnpt, zynpt, znnpt ! x,y components and norm of the vector: T point to North Pole
	141	REAL(wp) :: zxnpu, zynpu, znnpu ! x,y components and norm of the vector: U point to North Pole
	142	REAL(wp) :: zxnpv, zynpv, znnpv ! x,y components and norm of the vector: V point to North Pole
	143	REAL(wp) :: zxnpf, zynpf, znnpf ! x,y components and norm of the vector: F point to North Pole
	144	REAL(wp) :: zxvvt, zyvvt, znvvt ! x,y components and norm of the vector: between V points below and above a T point
	145	REAL(wp) :: zxffu, zyffu, znffu ! x,y components and norm of the vector: between F points below and above a U point
	146	REAL(wp) :: zxffv, zyffv, znffv ! x,y components and norm of the vector: between F points left and right a V point
	147	REAL(wp) :: zxuuf, zyuuf, znuuf ! x,y components and norm of the vector: between U points below and above a F point
[3]	148	!!----------------------------------------------------------------------
[6140]	149	!
[2715]	150	ALLOCATE( gsint(jpi,jpj), gcost(jpi,jpj), &
	151	& gsinu(jpi,jpj), gcosu(jpi,jpj), &
	152	& gsinv(jpi,jpj), gcosv(jpi,jpj), &
	153	& gsinf(jpi,jpj), gcosf(jpi,jpj), STAT=ierr )
[10425]	154	CALL mpp_sum( 'geo2ocean', ierr )
[6140]	155	IF( ierr /= 0 ) CALL ctl_stop( 'angle: unable to allocate arrays' )
	156	!
[3]	157	! ============================= !
	158	! Compute the cosinus and sinus !
	159	! ============================= !
	160	! (computation done on the north stereographic polar plane)
[6140]	161	!
[14215]	162	DO_2D( 0, 1, 0, 0 )
[12377]	163	!
	164	zlam = plamt(ji,jj) ! north pole direction & modulous (at t-point)
	165	zphi = pphit(ji,jj)
	166	zxnpt = 0. - 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	167	zynpt = 0. - 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	168	znnpt = zxnptzxnpt + zynptzynpt
	169	!
	170	zlam = plamu(ji,jj) ! north pole direction & modulous (at u-point)
	171	zphi = pphiu(ji,jj)
	172	zxnpu = 0. - 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	173	zynpu = 0. - 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	174	znnpu = zxnpuzxnpu + zynpuzynpu
	175	!
	176	zlam = plamv(ji,jj) ! north pole direction & modulous (at v-point)
	177	zphi = pphiv(ji,jj)
	178	zxnpv = 0. - 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	179	zynpv = 0. - 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	180	znnpv = zxnpvzxnpv + zynpvzynpv
	181	!
	182	zlam = plamf(ji,jj) ! north pole direction & modulous (at f-point)
	183	zphi = pphif(ji,jj)
	184	zxnpf = 0. - 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	185	zynpf = 0. - 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	186	znnpf = zxnpfzxnpf + zynpfzynpf
	187	!
	188	zlam = plamv(ji,jj ) ! j-direction: v-point segment direction (around t-point)
	189	zphi = pphiv(ji,jj )
	190	zlan = plamv(ji,jj-1)
	191	zphh = pphiv(ji,jj-1)
	192	zxvvt = 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	193	& - 2. * COS( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	194	zyvvt = 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	195	& - 2. * SIN( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	196	znvvt = SQRT( znnpt * ( zxvvtzxvvt + zyvvtzyvvt ) )
	197	znvvt = MAX( znvvt, 1.e-14 )
	198	!
	199	zlam = plamf(ji,jj ) ! j-direction: f-point segment direction (around u-point)
	200	zphi = pphif(ji,jj )
	201	zlan = plamf(ji,jj-1)
	202	zphh = pphif(ji,jj-1)
	203	zxffu = 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	204	& - 2. * COS( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	205	zyffu = 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	206	& - 2. * SIN( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	207	znffu = SQRT( znnpu * ( zxffuzxffu + zyffuzyffu ) )
	208	znffu = MAX( znffu, 1.e-14 )
	209	!
	210	zlam = plamf(ji ,jj) ! i-direction: f-point segment direction (around v-point)
	211	zphi = pphif(ji ,jj)
	212	zlan = plamf(ji-1,jj)
	213	zphh = pphif(ji-1,jj)
	214	zxffv = 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	215	& - 2. * COS( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	216	zyffv = 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	217	& - 2. * SIN( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	218	znffv = SQRT( znnpv * ( zxffvzxffv + zyffvzyffv ) )
	219	znffv = MAX( znffv, 1.e-14 )
	220	!
	221	zlam = plamu(ji,jj+1) ! j-direction: u-point segment direction (around f-point)
	222	zphi = pphiu(ji,jj+1)
	223	zlan = plamu(ji,jj )
	224	zphh = pphiu(ji,jj )
	225	zxuuf = 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	226	& - 2. * COS( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	227	zyuuf = 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	228	& - 2. * SIN( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	229	znuuf = SQRT( znnpf * ( zxuufzxuuf + zyuufzyuuf ) )
	230	znuuf = MAX( znuuf, 1.e-14 )
	231	!
	232	! ! cosinus and sinus using dot and cross products
	233	gsint(ji,jj) = ( zxnptzyvvt - zynptzxvvt ) / znvvt
	234	gcost(ji,jj) = ( zxnptzxvvt + zynptzyvvt ) / znvvt
	235	!
	236	gsinu(ji,jj) = ( zxnpuzyffu - zynpuzxffu ) / znffu
	237	gcosu(ji,jj) = ( zxnpuzxffu + zynpuzyffu ) / znffu
	238	!
	239	gsinf(ji,jj) = ( zxnpfzyuuf - zynpfzxuuf ) / znuuf
	240	gcosf(ji,jj) = ( zxnpfzxuuf + zynpfzyuuf ) / znuuf
	241	!
	242	gsinv(ji,jj) = ( zxnpvzxffv + zynpvzyffv ) / znffv
	243	gcosv(ji,jj) =-( zxnpvzyffv - zynpvzxffv ) / znffv ! (caution, rotation of 90 degres)
	244	!
	245	END_2D
[3]	246
	247	! =============== !
	248	! Geographic mesh !
	249	! =============== !
	250
[14215]	251	DO_2D( 0, 1, 0, 0 )
[12377]	252	IF( MOD( ABS( plamv(ji,jj) - plamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
	253	gsint(ji,jj) = 0.
	254	gcost(ji,jj) = 1.
	255	ENDIF
	256	IF( MOD( ABS( plamf(ji,jj) - plamf(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
	257	gsinu(ji,jj) = 0.
	258	gcosu(ji,jj) = 1.
	259	ENDIF
	260	IF( ABS( pphif(ji,jj) - pphif(ji-1,jj) ) < 1.e-8 ) THEN
	261	gsinv(ji,jj) = 0.
	262	gcosv(ji,jj) = 1.
	263	ENDIF
	264	IF( MOD( ABS( plamu(ji,jj) - plamu(ji,jj+1) ), 360. ) < 1.e-8 ) THEN
	265	gsinf(ji,jj) = 0.
	266	gcosf(ji,jj) = 1.
	267	ENDIF
	268	END_2D
[3]	269
	270	! =========================== !
	271	! Lateral boundary conditions !
	272	! =========================== !
[6140]	273	! ! lateral boundary cond.: T-, U-, V-, F-pts, sgn
[14433]	274	CALL lbc_lnk( 'geo2ocean', gcost, 'T', -1.0_wp, gsint, 'T', -1.0_wp, gcosu, 'U', -1.0_wp, gsinu, 'U', -1.0_wp, &
	275	& gcosv, 'V', -1.0_wp, gsinv, 'V', -1.0_wp, gcosf, 'F', -1.0_wp, gsinf, 'F', -1.0_wp )
[6140]	276	!
[3]	277	END SUBROUTINE angle
	278
	279
[6140]	280	SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid, pte, ptn )
[3]	281	!!----------------------------------------------------------------------
	282	!! * ROUTINE geo2oce *
	283	!!
	284	!! ** Purpose :
	285	!!
[6140]	286	!! ** Method : Change a vector from geocentric to east/north
[3]	287	!!
	288	!!----------------------------------------------------------------------
[1218]	289	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pxx, pyy, pzz
	290	CHARACTER(len=1) , INTENT(in ) :: cgrid
	291	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pte, ptn
[6140]	292	!
[2715]	293	REAL(wp), PARAMETER :: rpi = 3.141592653e0
[88]	294	REAL(wp), PARAMETER :: rad = rpi / 180.e0
[1218]	295	INTEGER :: ig !
[2715]	296	INTEGER :: ierr ! local integer
[1218]	297	!!----------------------------------------------------------------------
[6140]	298	!
[2715]	299	IF( .NOT. ALLOCATED( gsinlon ) ) THEN
	300	ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) , &
	301	& gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr )
[10425]	302	CALL mpp_sum( 'geo2ocean', ierr )
[4162]	303	IF( ierr /= 0 ) CALL ctl_stop('geo2oce: unable to allocate arrays' )
[2715]	304	ENDIF
[6140]	305	!
[1218]	306	SELECT CASE( cgrid)
[6140]	307	CASE ( 'T' )
	308	ig = 1
	309	IF( .NOT. linit(ig) ) THEN
	310	gsinlon(:,:,ig) = SIN( rad * glamt(:,:) )
	311	gcoslon(:,:,ig) = COS( rad * glamt(:,:) )
	312	gsinlat(:,:,ig) = SIN( rad * gphit(:,:) )
	313	gcoslat(:,:,ig) = COS( rad * gphit(:,:) )
	314	linit(ig) = .TRUE.
	315	ENDIF
	316	CASE ( 'U' )
	317	ig = 2
	318	IF( .NOT. linit(ig) ) THEN
	319	gsinlon(:,:,ig) = SIN( rad * glamu(:,:) )
	320	gcoslon(:,:,ig) = COS( rad * glamu(:,:) )
	321	gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) )
	322	gcoslat(:,:,ig) = COS( rad * gphiu(:,:) )
	323	linit(ig) = .TRUE.
	324	ENDIF
	325	CASE ( 'V' )
	326	ig = 3
	327	IF( .NOT. linit(ig) ) THEN
	328	gsinlon(:,:,ig) = SIN( rad * glamv(:,:) )
	329	gcoslon(:,:,ig) = COS( rad * glamv(:,:) )
	330	gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) )
	331	gcoslat(:,:,ig) = COS( rad * gphiv(:,:) )
	332	linit(ig) = .TRUE.
	333	ENDIF
	334	CASE ( 'F' )
	335	ig = 4
	336	IF( .NOT. linit(ig) ) THEN
	337	gsinlon(:,:,ig) = SIN( rad * glamf(:,:) )
	338	gcoslon(:,:,ig) = COS( rad * glamf(:,:) )
	339	gsinlat(:,:,ig) = SIN( rad * gphif(:,:) )
	340	gcoslat(:,:,ig) = COS( rad * gphif(:,:) )
	341	linit(ig) = .TRUE.
	342	ENDIF
	343	CASE default
	344	WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid
	345	CALL ctl_stop( ctmp1 )
[1218]	346	END SELECT
[6140]	347	!
[2715]	348	pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy
	349	ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx &
[6140]	350	& - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy &
	351	& + gcoslat(:,:,ig) * pzz
[1218]	352	!
	353	END SUBROUTINE geo2oce
	354
[6140]	355
	356	SUBROUTINE oce2geo ( pte, ptn, cgrid, pxx , pyy , pzz )
[1218]	357	!!----------------------------------------------------------------------
	358	!! * ROUTINE oce2geo *
	359	!!
	360	!! ** Purpose :
	361	!!
	362	!! ** Method : Change vector from east/north to geocentric
	363	!!
[6140]	364	!! History : ! (A. Caubel) oce2geo - Original code
[1218]	365	!!----------------------------------------------------------------------
	366	REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) :: pte, ptn
	367	CHARACTER(len=1) , INTENT( IN ) :: cgrid
	368	REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ) :: pxx , pyy , pzz
	369	!!
	370	REAL(wp), PARAMETER :: rpi = 3.141592653E0
	371	REAL(wp), PARAMETER :: rad = rpi / 180.e0
[3]	372	INTEGER :: ig !
[2715]	373	INTEGER :: ierr ! local integer
[3]	374	!!----------------------------------------------------------------------
	375
[4162]	376	IF( .NOT. ALLOCATED( gsinlon ) ) THEN
[2715]	377	ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) , &
	378	& gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr )
[10425]	379	CALL mpp_sum( 'geo2ocean', ierr )
[4162]	380	IF( ierr /= 0 ) CALL ctl_stop('oce2geo: unable to allocate arrays' )
[2715]	381	ENDIF
	382
[3]	383	SELECT CASE( cgrid)
[1226]	384	CASE ( 'T' )
	385	ig = 1
	386	IF( .NOT. linit(ig) ) THEN
[2715]	387	gsinlon(:,:,ig) = SIN( rad * glamt(:,:) )
	388	gcoslon(:,:,ig) = COS( rad * glamt(:,:) )
	389	gsinlat(:,:,ig) = SIN( rad * gphit(:,:) )
	390	gcoslat(:,:,ig) = COS( rad * gphit(:,:) )
[1226]	391	linit(ig) = .TRUE.
	392	ENDIF
	393	CASE ( 'U' )
	394	ig = 2
	395	IF( .NOT. linit(ig) ) THEN
[2715]	396	gsinlon(:,:,ig) = SIN( rad * glamu(:,:) )
	397	gcoslon(:,:,ig) = COS( rad * glamu(:,:) )
	398	gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) )
	399	gcoslat(:,:,ig) = COS( rad * gphiu(:,:) )
[1226]	400	linit(ig) = .TRUE.
	401	ENDIF
	402	CASE ( 'V' )
	403	ig = 3
	404	IF( .NOT. linit(ig) ) THEN
[2715]	405	gsinlon(:,:,ig) = SIN( rad * glamv(:,:) )
	406	gcoslon(:,:,ig) = COS( rad * glamv(:,:) )
	407	gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) )
	408	gcoslat(:,:,ig) = COS( rad * gphiv(:,:) )
[1226]	409	linit(ig) = .TRUE.
	410	ENDIF
	411	CASE ( 'F' )
	412	ig = 4
	413	IF( .NOT. linit(ig) ) THEN
[2715]	414	gsinlon(:,:,ig) = SIN( rad * glamf(:,:) )
	415	gcoslon(:,:,ig) = COS( rad * glamf(:,:) )
	416	gsinlat(:,:,ig) = SIN( rad * gphif(:,:) )
	417	gcoslat(:,:,ig) = COS( rad * gphif(:,:) )
[1226]	418	linit(ig) = .TRUE.
	419	ENDIF
	420	CASE default
	421	WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid
[474]	422	CALL ctl_stop( ctmp1 )
[1226]	423	END SELECT
[6140]	424	!
	425	pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn
	426	pyy = gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn
	427	pzz = gcoslat(:,:,ig) * ptn
	428	!
[1218]	429	END SUBROUTINE oce2geo
[3]	430
	431
[6140]	432	SUBROUTINE obs_rot( psinu, pcosu, psinv, pcosv )
[3]	433	!!----------------------------------------------------------------------
[2528]	434	!! * ROUTINE obs_rot *
	435	!!
	436	!! ** Purpose : Copy gsinu, gcosu, gsinv and gsinv
	437	!! to input data for rotations of
	438	!! current at observation points
	439	!!
[6140]	440	!! History : 9.2 ! 09-02 (K. Mogensen)
[2528]	441	!!----------------------------------------------------------------------
[2715]	442	REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ):: psinu, pcosu, psinv, pcosv ! copy of data
[2528]	443	!!----------------------------------------------------------------------
[6140]	444	!
[2528]	445	! Initialization of gsin* and gcos* at first call
	446	! -----------------------------------------------
	447	IF( lmust_init ) THEN
	448	IF(lwp) WRITE(numout,*)
	449	IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched'
	450	IF(lwp) WRITE(numout,*) ' ~~~~~~~ coordinate transformation'
[10370]	451	CALL angle( glamt, gphit, glamu, gphiu, glamv, gphiv, glamf, gphif ) ! initialization of the transformation
[2528]	452	lmust_init = .FALSE.
	453	ENDIF
[6140]	454	!
[2528]	455	psinu(:,:) = gsinu(:,:)
	456	pcosu(:,:) = gcosu(:,:)
	457	psinv(:,:) = gsinv(:,:)
	458	pcosv(:,:) = gcosv(:,:)
[6140]	459	!
[2528]	460	END SUBROUTINE obs_rot
	461
[3]	462	!!======================================================================
	463	END MODULE geo2ocean

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

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

Context Navigation

source: NEMO/trunk/src/OCE/SBC/geo2ocean.F90

Download in other formats: