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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90 @ 7702

Last change on this file since 7702 was 7698, checked in by mocavero, 7 years ago
update trunk with OpenMP parallelization
Property svn:keywords set to `Id`
File size: 25.0 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
[7698]	68	INTEGER :: ji, jj ! dummy loop indices
[6140]	69	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: prot
[672]	70	!!----------------------------------------------------------------------
[6140]	71	!
	72	IF( lmust_init ) THEN ! at 1st call only: set gsin. & gcos.
[3]	73	IF(lwp) WRITE(numout,*)
[6140]	74	IF(lwp) WRITE(numout,*) ' rot_rep: coordinate transformation : geographic <==> model (i,j)-components'
	75	IF(lwp) WRITE(numout,*) ' ~~~~~~~~ '
[2715]	76	!
[3]	77	CALL angle ! initialization of the transformation
[672]	78	lmust_init = .FALSE.
[3]	79	ENDIF
[6140]	80	!
	81	SELECT CASE( cdtodo ) ! type of rotation
	82	!
	83	CASE( 'en->i' ) ! east-north to i-component
[672]	84	SELECT CASE (cd_type)
[7698]	85	CASE ('T')
	86	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	87	DO jj = 1, jpj
	88	DO ji = 1, jpi
	89	prot(ji,jj) = pxin(ji,jj) * gcost(ji,jj) + pyin(ji,jj) * gsint(ji,jj)
	90	END DO
	91	END DO
	92	CASE ('U')
	93	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	94	DO jj = 1, jpj
	95	DO ji = 1, jpi
	96	prot(ji,jj) = pxin(ji,jj) * gcosu(ji,jj) + pyin(ji,jj) * gsinu(ji,jj)
	97	END DO
	98	END DO
	99	CASE ('V')
	100	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	101	DO jj = 1, jpj
	102	DO ji = 1, jpi
	103	prot(ji,jj) = pxin(ji,jj) * gcosv(ji,jj) + pyin(ji,jj) * gsinv(ji,jj)
	104	END DO
	105	END DO
	106	CASE ('F')
	107	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	108	DO jj = 1, jpj
	109	DO ji = 1, jpi
	110	prot(ji,jj) = pxin(ji,jj) * gcosf(ji,jj) + pyin(ji,jj) * gsinf(ji,jj)
	111	END DO
	112	END DO
[672]	113	CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
	114	END SELECT
[6140]	115	CASE ('en->j') ! east-north to j-component
[672]	116	SELECT CASE (cd_type)
[7698]	117	CASE ('T')
	118	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	119	DO jj = 1, jpj
	120	DO ji = 1, jpi
	121	prot(ji,jj) = pyin(ji,jj) * gcost(ji,jj) - pxin(ji,jj) * gsint(ji,jj)
	122	END DO
	123	END DO
	124	CASE ('U')
	125	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	126	DO jj = 1, jpj
	127	DO ji = 1, jpi
	128	prot(ji,jj) = pyin(ji,jj) * gcosu(ji,jj) - pxin(ji,jj) * gsinu(ji,jj)
	129	END DO
	130	END DO
	131	CASE ('V')
	132	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	133	DO jj = 1, jpj
	134	DO ji = 1, jpi
	135	prot(ji,jj) = pyin(ji,jj) * gcosv(ji,jj) - pxin(ji,jj) * gsinv(ji,jj)
	136	END DO
	137	END DO
	138	CASE ('F')
	139	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	140	DO jj = 1, jpj
	141	DO ji = 1, jpi
	142	prot(ji,jj) = pyin(ji,jj) * gcosf(ji,jj) - pxin(ji,jj) * gsinf(ji,jj)
	143	END DO
	144	END DO
[672]	145	CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
	146	END SELECT
[6140]	147	CASE ('ij->e') ! (i,j)-components to east
[672]	148	SELECT CASE (cd_type)
[7698]	149	CASE ('T')
	150	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	151	DO jj = 1, jpj
	152	DO ji = 1, jpi
	153	prot(ji,jj) = pxin(ji,jj) * gcost(ji,jj) - pyin(ji,jj) * gsint(ji,jj)
	154	END DO
	155	END DO
	156	CASE ('U')
	157	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	158	DO jj = 1, jpj
	159	DO ji = 1, jpi
	160	prot(ji,jj) = pxin(ji,jj) * gcosu(ji,jj) - pyin(ji,jj) * gsinu(ji,jj)
	161	END DO
	162	END DO
	163	CASE ('V')
	164	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	165	DO jj = 1, jpj
	166	DO ji = 1, jpi
	167	prot(ji,jj) = pxin(ji,jj) * gcosv(ji,jj) - pyin(ji,jj) * gsinv(ji,jj)
	168	END DO
	169	END DO
	170	CASE ('F')
	171	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	172	DO jj = 1, jpj
	173	DO ji = 1, jpi
	174	prot(ji,jj) = pxin(ji,jj) * gcosf(ji,jj) - pyin(ji,jj) * gsinf(ji,jj)
	175	END DO
	176	END DO
[672]	177	CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
	178	END SELECT
[6140]	179	CASE ('ij->n') ! (i,j)-components to north
[672]	180	SELECT CASE (cd_type)
[7698]	181	CASE ('T')
	182	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	183	DO jj = 1, jpj
	184	DO ji = 1, jpi
	185	prot(ji,jj) = pyin(ji,jj) * gcost(ji,jj) + pxin(ji,jj) * gsint(ji,jj)
	186	END DO
	187	END DO
	188	CASE ('U')
	189	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	190	DO jj = 1, jpj
	191	DO ji = 1, jpi
	192	prot(ji,jj) = pyin(ji,jj) * gcosu(ji,jj) + pxin(ji,jj) * gsinu(ji,jj)
	193	END DO
	194	END DO
	195	CASE ('V')
	196	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	197	DO jj = 1, jpj
	198	DO ji = 1, jpi
	199	prot(ji,jj) = pyin(ji,jj) * gcosv(ji,jj) + pxin(ji,jj) * gsinv(ji,jj)
	200	END DO
	201	END DO
	202	CASE ('F')
	203	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	204	DO jj = 1, jpj
	205	DO ji = 1, jpi
	206	prot(ji,jj) = pyin(ji,jj) * gcosf(ji,jj) + pxin(ji,jj) * gsinf(ji,jj)
	207	END DO
	208	END DO
[672]	209	CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
	210	END SELECT
	211	CASE DEFAULT ; CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' )
[6140]	212	!
[672]	213	END SELECT
[6140]	214	!
[685]	215	END SUBROUTINE rot_rep
[3]	216
	217
	218	SUBROUTINE angle
	219	!!----------------------------------------------------------------------
	220	!! * ROUTINE angle *
	221	!!
[672]	222	!! ** Purpose : Compute angles between model grid lines and the North direction
[3]	223	!!
[6140]	224	!! ** Method : sinus and cosinus of the angle between the north-south axe
	225	!! and the j-direction at t, u, v and f-points
	226	!! dot and cross products are used to obtain cos and sin, resp.
[3]	227	!!
[6140]	228	!! ** Action : - gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf
[3]	229	!!----------------------------------------------------------------------
[6140]	230	INTEGER :: ji, jj ! dummy loop indices
	231	INTEGER :: ierr ! local integer
	232	REAL(wp) :: zlam, zphi ! local scalars
	233	REAL(wp) :: zlan, zphh ! - -
	234	REAL(wp) :: zxnpt, zynpt, znnpt ! x,y components and norm of the vector: T point to North Pole
	235	REAL(wp) :: zxnpu, zynpu, znnpu ! x,y components and norm of the vector: U point to North Pole
	236	REAL(wp) :: zxnpv, zynpv, znnpv ! x,y components and norm of the vector: V point to North Pole
	237	REAL(wp) :: zxnpf, zynpf, znnpf ! x,y components and norm of the vector: F point to North Pole
	238	REAL(wp) :: zxvvt, zyvvt, znvvt ! x,y components and norm of the vector: between V points below and above a T point
	239	REAL(wp) :: zxffu, zyffu, znffu ! x,y components and norm of the vector: between F points below and above a U point
	240	REAL(wp) :: zxffv, zyffv, znffv ! x,y components and norm of the vector: between F points left and right a V point
	241	REAL(wp) :: zxuuf, zyuuf, znuuf ! x,y components and norm of the vector: between U points below and above a F point
[3]	242	!!----------------------------------------------------------------------
[6140]	243	!
[2715]	244	ALLOCATE( gsint(jpi,jpj), gcost(jpi,jpj), &
	245	& gsinu(jpi,jpj), gcosu(jpi,jpj), &
	246	& gsinv(jpi,jpj), gcosv(jpi,jpj), &
	247	& gsinf(jpi,jpj), gcosf(jpi,jpj), STAT=ierr )
	248	IF(lk_mpp) CALL mpp_sum( ierr )
[6140]	249	IF( ierr /= 0 ) CALL ctl_stop( 'angle: unable to allocate arrays' )
	250	!
[3]	251	! ============================= !
	252	! Compute the cosinus and sinus !
	253	! ============================= !
	254	! (computation done on the north stereographic polar plane)
[6140]	255	!
[7698]	256	!$OMP PARALLEL
	257	!$OMP DO schedule(static) private(jj,ji,zlam,zphi,zxnpt,zynpt,znnpt,zxnpu,zynpu,znnpu,zxnpv,zynpv,znnpv,zxnpf) &
	258	!$OMP& private(zynpf,znnpf,zlan,zphh,zxvvt,zyvvt,znvvt,zxffu,zyffu,znffu,zxffv,zyffv,znffv,zxuuf,zyuuf,znuuf)
[672]	259	DO jj = 2, jpjm1
[3]	260	DO ji = fs_2, jpi ! vector opt.
[6140]	261	!
	262	zlam = glamt(ji,jj) ! north pole direction & modulous (at t-point)
[672]	263	zphi = gphit(ji,jj)
	264	zxnpt = 0. - 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	265	zynpt = 0. - 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	266	znnpt = zxnptzxnpt + zynptzynpt
[6140]	267	!
	268	zlam = glamu(ji,jj) ! north pole direction & modulous (at u-point)
[3]	269	zphi = gphiu(ji,jj)
	270	zxnpu = 0. - 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	271	zynpu = 0. - 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	272	znnpu = zxnpuzxnpu + zynpuzynpu
[6140]	273	!
	274	zlam = glamv(ji,jj) ! north pole direction & modulous (at v-point)
[3]	275	zphi = gphiv(ji,jj)
	276	zxnpv = 0. - 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	277	zynpv = 0. - 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	278	znnpv = zxnpvzxnpv + zynpvzynpv
[6140]	279	!
	280	zlam = glamf(ji,jj) ! north pole direction & modulous (at f-point)
[672]	281	zphi = gphif(ji,jj)
	282	zxnpf = 0. - 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	283	zynpf = 0. - 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	284	znnpf = zxnpfzxnpf + zynpfzynpf
[6140]	285	!
	286	zlam = glamv(ji,jj ) ! j-direction: v-point segment direction (around t-point)
[672]	287	zphi = gphiv(ji,jj )
	288	zlan = glamv(ji,jj-1)
	289	zphh = gphiv(ji,jj-1)
	290	zxvvt = 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	291	& - 2. * COS( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	292	zyvvt = 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	293	& - 2. * SIN( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	294	znvvt = SQRT( znnpt * ( zxvvtzxvvt + zyvvtzyvvt ) )
	295	znvvt = MAX( znvvt, 1.e-14 )
[6140]	296	!
	297	zlam = glamf(ji,jj ) ! j-direction: f-point segment direction (around u-point)
[3]	298	zphi = gphif(ji,jj )
	299	zlan = glamf(ji,jj-1)
	300	zphh = gphif(ji,jj-1)
	301	zxffu = 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	302	& - 2. * COS( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	303	zyffu = 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	304	& - 2. * SIN( radzlan ) TAN( rpi/4. - rad*zphh/2. )
[672]	305	znffu = SQRT( znnpu * ( zxffuzxffu + zyffuzyffu ) )
	306	znffu = MAX( znffu, 1.e-14 )
[6140]	307	!
	308	zlam = glamf(ji ,jj) ! i-direction: f-point segment direction (around v-point)
[3]	309	zphi = gphif(ji ,jj)
	310	zlan = glamf(ji-1,jj)
	311	zphh = gphif(ji-1,jj)
	312	zxffv = 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	313	& - 2. * COS( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	314	zyffv = 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	315	& - 2. * SIN( radzlan ) TAN( rpi/4. - rad*zphh/2. )
[672]	316	znffv = SQRT( znnpv * ( zxffvzxffv + zyffvzyffv ) )
	317	znffv = MAX( znffv, 1.e-14 )
[6140]	318	!
	319	zlam = glamu(ji,jj+1) ! j-direction: u-point segment direction (around f-point)
[672]	320	zphi = gphiu(ji,jj+1)
	321	zlan = glamu(ji,jj )
	322	zphh = gphiu(ji,jj )
	323	zxuuf = 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	324	& - 2. * COS( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	325	zyuuf = 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	326	& - 2. * SIN( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	327	znuuf = SQRT( znnpf * ( zxuufzxuuf + zyuufzyuuf ) )
	328	znuuf = MAX( znuuf, 1.e-14 )
[6140]	329	!
	330	! ! cosinus and sinus using dot and cross products
[672]	331	gsint(ji,jj) = ( zxnptzyvvt - zynptzxvvt ) / znvvt
	332	gcost(ji,jj) = ( zxnptzxvvt + zynptzyvvt ) / znvvt
[6140]	333	!
[672]	334	gsinu(ji,jj) = ( zxnpuzyffu - zynpuzxffu ) / znffu
	335	gcosu(ji,jj) = ( zxnpuzxffu + zynpuzyffu ) / znffu
[6140]	336	!
[672]	337	gsinf(ji,jj) = ( zxnpfzyuuf - zynpfzxuuf ) / znuuf
	338	gcosf(ji,jj) = ( zxnpfzxuuf + zynpfzyuuf ) / znuuf
[6140]	339	!
[672]	340	gsinv(ji,jj) = ( zxnpvzxffv + zynpvzyffv ) / znffv
[6140]	341	gcosv(ji,jj) =-( zxnpvzyffv - zynpvzxffv ) / znffv ! (caution, rotation of 90 degres)
	342	!
[3]	343	END DO
	344	END DO
	345
	346	! =============== !
	347	! Geographic mesh !
	348	! =============== !
	349
[7698]	350	!$OMP DO schedule(static) private(jj,ji)
[672]	351	DO jj = 2, jpjm1
[3]	352	DO ji = fs_2, jpi ! vector opt.
[672]	353	IF( MOD( ABS( glamv(ji,jj) - glamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
	354	gsint(ji,jj) = 0.
	355	gcost(ji,jj) = 1.
	356	ENDIF
	357	IF( MOD( ABS( glamf(ji,jj) - glamf(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
[3]	358	gsinu(ji,jj) = 0.
	359	gcosu(ji,jj) = 1.
	360	ENDIF
[672]	361	IF( ABS( gphif(ji,jj) - gphif(ji-1,jj) ) < 1.e-8 ) THEN
[3]	362	gsinv(ji,jj) = 0.
	363	gcosv(ji,jj) = 1.
	364	ENDIF
[672]	365	IF( MOD( ABS( glamu(ji,jj) - glamu(ji,jj+1) ), 360. ) < 1.e-8 ) THEN
	366	gsinf(ji,jj) = 0.
	367	gcosf(ji,jj) = 1.
	368	ENDIF
[3]	369	END DO
	370	END DO
[7698]	371	!$OMP END DO NOWAIT
	372	!$OMP END PARALLEL
[3]	373
	374	! =========================== !
	375	! Lateral boundary conditions !
	376	! =========================== !
[6140]	377	! ! lateral boundary cond.: T-, U-, V-, F-pts, sgn
[1833]	378	CALL lbc_lnk( gcost, 'T', -1. ) ; CALL lbc_lnk( gsint, 'T', -1. )
	379	CALL lbc_lnk( gcosu, 'U', -1. ) ; CALL lbc_lnk( gsinu, 'U', -1. )
	380	CALL lbc_lnk( gcosv, 'V', -1. ) ; CALL lbc_lnk( gsinv, 'V', -1. )
	381	CALL lbc_lnk( gcosf, 'F', -1. ) ; CALL lbc_lnk( gsinf, 'F', -1. )
[6140]	382	!
[3]	383	END SUBROUTINE angle
	384
	385
[6140]	386	SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid, pte, ptn )
[3]	387	!!----------------------------------------------------------------------
	388	!! * ROUTINE geo2oce *
	389	!!
	390	!! ** Purpose :
	391	!!
[6140]	392	!! ** Method : Change a vector from geocentric to east/north
[3]	393	!!
	394	!!----------------------------------------------------------------------
[1218]	395	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pxx, pyy, pzz
	396	CHARACTER(len=1) , INTENT(in ) :: cgrid
	397	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pte, ptn
[6140]	398	!
[2715]	399	REAL(wp), PARAMETER :: rpi = 3.141592653e0
[88]	400	REAL(wp), PARAMETER :: rad = rpi / 180.e0
[1218]	401	INTEGER :: ig !
[2715]	402	INTEGER :: ierr ! local integer
[1218]	403	!!----------------------------------------------------------------------
[6140]	404	!
[2715]	405	IF( .NOT. ALLOCATED( gsinlon ) ) THEN
	406	ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) , &
	407	& gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr )
	408	IF( lk_mpp ) CALL mpp_sum( ierr )
[4162]	409	IF( ierr /= 0 ) CALL ctl_stop('geo2oce: unable to allocate arrays' )
[2715]	410	ENDIF
[6140]	411	!
[1218]	412	SELECT CASE( cgrid)
[6140]	413	CASE ( 'T' )
	414	ig = 1
	415	IF( .NOT. linit(ig) ) THEN
	416	gsinlon(:,:,ig) = SIN( rad * glamt(:,:) )
	417	gcoslon(:,:,ig) = COS( rad * glamt(:,:) )
	418	gsinlat(:,:,ig) = SIN( rad * gphit(:,:) )
	419	gcoslat(:,:,ig) = COS( rad * gphit(:,:) )
	420	linit(ig) = .TRUE.
	421	ENDIF
	422	CASE ( 'U' )
	423	ig = 2
	424	IF( .NOT. linit(ig) ) THEN
	425	gsinlon(:,:,ig) = SIN( rad * glamu(:,:) )
	426	gcoslon(:,:,ig) = COS( rad * glamu(:,:) )
	427	gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) )
	428	gcoslat(:,:,ig) = COS( rad * gphiu(:,:) )
	429	linit(ig) = .TRUE.
	430	ENDIF
	431	CASE ( 'V' )
	432	ig = 3
	433	IF( .NOT. linit(ig) ) THEN
	434	gsinlon(:,:,ig) = SIN( rad * glamv(:,:) )
	435	gcoslon(:,:,ig) = COS( rad * glamv(:,:) )
	436	gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) )
	437	gcoslat(:,:,ig) = COS( rad * gphiv(:,:) )
	438	linit(ig) = .TRUE.
	439	ENDIF
	440	CASE ( 'F' )
	441	ig = 4
	442	IF( .NOT. linit(ig) ) THEN
	443	gsinlon(:,:,ig) = SIN( rad * glamf(:,:) )
	444	gcoslon(:,:,ig) = COS( rad * glamf(:,:) )
	445	gsinlat(:,:,ig) = SIN( rad * gphif(:,:) )
	446	gcoslat(:,:,ig) = COS( rad * gphif(:,:) )
	447	linit(ig) = .TRUE.
	448	ENDIF
	449	CASE default
	450	WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid
	451	CALL ctl_stop( ctmp1 )
[1218]	452	END SELECT
[6140]	453	!
[2715]	454	pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy
	455	ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx &
[6140]	456	& - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy &
	457	& + gcoslat(:,:,ig) * pzz
[1218]	458	!
	459	END SUBROUTINE geo2oce
	460
[6140]	461
	462	SUBROUTINE oce2geo ( pte, ptn, cgrid, pxx , pyy , pzz )
[1218]	463	!!----------------------------------------------------------------------
	464	!! * ROUTINE oce2geo *
	465	!!
	466	!! ** Purpose :
	467	!!
	468	!! ** Method : Change vector from east/north to geocentric
	469	!!
[6140]	470	!! History : ! (A. Caubel) oce2geo - Original code
[1218]	471	!!----------------------------------------------------------------------
	472	REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) :: pte, ptn
	473	CHARACTER(len=1) , INTENT( IN ) :: cgrid
	474	REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ) :: pxx , pyy , pzz
	475	!!
	476	REAL(wp), PARAMETER :: rpi = 3.141592653E0
	477	REAL(wp), PARAMETER :: rad = rpi / 180.e0
[3]	478	INTEGER :: ig !
[2715]	479	INTEGER :: ierr ! local integer
[3]	480	!!----------------------------------------------------------------------
	481
[4162]	482	IF( .NOT. ALLOCATED( gsinlon ) ) THEN
[2715]	483	ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) , &
	484	& gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr )
	485	IF( lk_mpp ) CALL mpp_sum( ierr )
[4162]	486	IF( ierr /= 0 ) CALL ctl_stop('oce2geo: unable to allocate arrays' )
[2715]	487	ENDIF
	488
[3]	489	SELECT CASE( cgrid)
[1226]	490	CASE ( 'T' )
	491	ig = 1
	492	IF( .NOT. linit(ig) ) THEN
[2715]	493	gsinlon(:,:,ig) = SIN( rad * glamt(:,:) )
	494	gcoslon(:,:,ig) = COS( rad * glamt(:,:) )
	495	gsinlat(:,:,ig) = SIN( rad * gphit(:,:) )
	496	gcoslat(:,:,ig) = COS( rad * gphit(:,:) )
[1226]	497	linit(ig) = .TRUE.
	498	ENDIF
	499	CASE ( 'U' )
	500	ig = 2
	501	IF( .NOT. linit(ig) ) THEN
[2715]	502	gsinlon(:,:,ig) = SIN( rad * glamu(:,:) )
	503	gcoslon(:,:,ig) = COS( rad * glamu(:,:) )
	504	gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) )
	505	gcoslat(:,:,ig) = COS( rad * gphiu(:,:) )
[1226]	506	linit(ig) = .TRUE.
	507	ENDIF
	508	CASE ( 'V' )
	509	ig = 3
	510	IF( .NOT. linit(ig) ) THEN
[2715]	511	gsinlon(:,:,ig) = SIN( rad * glamv(:,:) )
	512	gcoslon(:,:,ig) = COS( rad * glamv(:,:) )
	513	gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) )
	514	gcoslat(:,:,ig) = COS( rad * gphiv(:,:) )
[1226]	515	linit(ig) = .TRUE.
	516	ENDIF
	517	CASE ( 'F' )
	518	ig = 4
	519	IF( .NOT. linit(ig) ) THEN
[2715]	520	gsinlon(:,:,ig) = SIN( rad * glamf(:,:) )
	521	gcoslon(:,:,ig) = COS( rad * glamf(:,:) )
	522	gsinlat(:,:,ig) = SIN( rad * gphif(:,:) )
	523	gcoslat(:,:,ig) = COS( rad * gphif(:,:) )
[1226]	524	linit(ig) = .TRUE.
	525	ENDIF
	526	CASE default
	527	WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid
[474]	528	CALL ctl_stop( ctmp1 )
[1226]	529	END SELECT
[6140]	530	!
	531	pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn
	532	pyy = gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn
	533	pzz = gcoslat(:,:,ig) * ptn
	534	!
[1218]	535	END SUBROUTINE oce2geo
[3]	536
	537
[6140]	538	SUBROUTINE obs_rot( psinu, pcosu, psinv, pcosv )
[3]	539	!!----------------------------------------------------------------------
[2528]	540	!! * ROUTINE obs_rot *
	541	!!
	542	!! ** Purpose : Copy gsinu, gcosu, gsinv and gsinv
	543	!! to input data for rotations of
	544	!! current at observation points
	545	!!
[6140]	546	!! History : 9.2 ! 09-02 (K. Mogensen)
[2528]	547	!!----------------------------------------------------------------------
[2715]	548	REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ):: psinu, pcosu, psinv, pcosv ! copy of data
[2528]	549	!!----------------------------------------------------------------------
[6140]	550	!
[2528]	551	! Initialization of gsin* and gcos* at first call
	552	! -----------------------------------------------
	553	IF( lmust_init ) THEN
	554	IF(lwp) WRITE(numout,*)
	555	IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched'
	556	IF(lwp) WRITE(numout,*) ' ~~~~~~~ coordinate transformation'
	557	CALL angle ! initialization of the transformation
	558	lmust_init = .FALSE.
	559	ENDIF
[6140]	560	!
[2528]	561	psinu(:,:) = gsinu(:,:)
	562	pcosu(:,:) = gcosu(:,:)
	563	psinv(:,:) = gsinv(:,:)
	564	pcosv(:,:) = gcosv(:,:)
[6140]	565	!
[2528]	566	END SUBROUTINE obs_rot
	567
[3]	568	!!======================================================================
	569	END MODULE geo2ocean

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

Download in other formats: