Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

geo2ocean.F90 @ 11105

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

Download in other formats: