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.

geo2ocean.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90 @ 9213

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

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

Download in other formats: