Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

geo2ocean.F90 @ 3294

Last change on this file since 3294 was 1884, checked in by rblod, 14 years ago
Light adaptation of NEMO direct model routine to handle TAM
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 25.2 KB

Rev	Line
[3]	1	MODULE geo2ocean
	2	!!======================================================================
	3	!! * MODULE geo2ocean *
[144]	4	!! Ocean mesh : ???
[1884]	5	!!======================================================================
	6	!! History : OPA ! 07-1996 (O. Marti) Original code
	7	!! NEMO 1.0 ! 02-2008 (G. Madec) F90: Free form
	8	!! 3.0 !
	9	!!----------------------------------------------------------------------
[3]	10
	11	!!----------------------------------------------------------------------
	12	!! repcmo :
	13	!! angle :
	14	!! geo2oce :
	15	!! repere : old routine suppress it ???
	16	!!----------------------------------------------------------------------
	17	USE dom_oce ! mesh and scale factors
	18	USE phycst ! physical constants
	19	USE in_out_manager ! I/O manager
	20	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	21
	22	IMPLICIT NONE
[1884]	23	PRIVATE
[3]	24
[1884]	25	PUBLIC rot_rep, repcmo, repere, geo2oce, oce2geo ! only rot_rep should be used
[672]	26	! repcmo and repere are keep only for compatibility.
	27	! they are only a useless overlay of rot_rep
[1884]	28	PUBLIC obs_rot
[3]	29
	30	REAL(wp), DIMENSION(jpi,jpj) :: &
[672]	31	gsint, gcost, & ! cos/sin between model grid lines and NP direction at T point
	32	gsinu, gcosu, & ! cos/sin between model grid lines and NP direction at U point
	33	gsinv, gcosv, & ! cos/sin between model grid lines and NP direction at V point
	34	gsinf, gcosf ! cos/sin between model grid lines and NP direction at F point
[3]	35
[672]	36	LOGICAL :: lmust_init = .TRUE. !: used to initialize the cos/sin variables (se above)
	37
[1884]	38	!! * Substitutions
[3]	39	# include "vectopt_loop_substitute.h90"
[1884]	40	!!----------------------------------------------------------------------
	41	!! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
	42	!! $Id$
	43	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
	44	!!----------------------------------------------------------------------
[3]	45
	46	CONTAINS
	47
	48	SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1, &
[685]	49	px2 , py2 )
[3]	50	!!----------------------------------------------------------------------
	51	!! * ROUTINE repcmo *
	52	!!
	53	!! ** Purpose : Change vector componantes from a geographic grid to a
	54	!! stretched coordinates grid.
	55	!!
	56	!! ** Method : Initialization of arrays at the first call.
	57	!!
[1884]	58	!! ** Action : - px2 : first componante (defined at u point)
[3]	59	!! - py2 : second componante (defined at v point)
	60	!!----------------------------------------------------------------------
[1884]	61	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pxu1, pyu1 ! geographic vector componantes at u-point
	62	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pxv1, pyv1 ! geographic vector componantes at v-point
	63	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: px2 ! i-componante (defined at u-point)
	64	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: py2 ! j-componante (defined at v-point)
[3]	65	!!----------------------------------------------------------------------
[672]	66
	67	! Change from geographic to stretched coordinate
	68	! ----------------------------------------------
[685]	69	CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 )
	70	CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 )
[672]	71
	72	END SUBROUTINE repcmo
[3]	73
	74
[685]	75	SUBROUTINE rot_rep ( pxin, pyin, cd_type, cdtodo, prot )
[672]	76	!!----------------------------------------------------------------------
	77	!! * ROUTINE rot_rep *
	78	!!
	79	!! ** Purpose : Rotate the Repere: Change vector componantes between
	80	!! geographic grid <--> stretched coordinates grid.
	81	!!
	82	!! History :
	83	!! 9.2 ! 07-04 (S. Masson)
	84	!! (O. Marti ) Original code (repere and repcmo)
	85	!!----------------------------------------------------------------------
	86	REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) :: pxin, pyin ! vector componantes
	87	CHARACTER(len=1), INTENT( IN ) :: cd_type ! define the nature of pt2d array grid-points
	88	CHARACTER(len=5), INTENT( IN ) :: cdtodo ! specify the work to do:
	89	!! ! 'en->i' east-north componantes to model i componante
	90	!! ! 'en->j' east-north componantes to model j componante
	91	!! ! 'ij->e' model i-j componantes to east componante
	92	!! ! 'ij->n' model i-j componantes to east componante
[685]	93	REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: prot
[672]	94
	95	!!----------------------------------------------------------------------
	96
[3]	97	! Initialization of gsin* and gcos* at first call
	98	! -----------------------------------------------
	99
[672]	100	IF( lmust_init ) THEN
[3]	101	IF(lwp) WRITE(numout,*)
[672]	102	IF(lwp) WRITE(numout,*) ' rot_rep : geographic <--> stretched'
	103	IF(lwp) WRITE(numout,*) ' ~~~~~ coordinate transformation'
[3]	104
	105	CALL angle ! initialization of the transformation
[672]	106	lmust_init = .FALSE.
	107
[3]	108	ENDIF
	109
[672]	110	SELECT CASE (cdtodo)
	111	CASE ('en->i') ! 'en->i' est-north componantes to model i componante
	112	SELECT CASE (cd_type)
[685]	113	CASE ('T') ; prot(:,:) = pxin(:,:) * gcost(:,:) + pyin(:,:) * gsint(:,:)
	114	CASE ('U') ; prot(:,:) = pxin(:,:) * gcosu(:,:) + pyin(:,:) * gsinu(:,:)
	115	CASE ('V') ; prot(:,:) = pxin(:,:) * gcosv(:,:) + pyin(:,:) * gsinv(:,:)
	116	CASE ('F') ; prot(:,:) = pxin(:,:) * gcosf(:,:) + pyin(:,:) * gsinf(:,:)
[672]	117	CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
	118	END SELECT
	119	CASE ('en->j') ! 'en->j' est-north componantes to model j componante
	120	SELECT CASE (cd_type)
[685]	121	CASE ('T') ; prot(:,:) = pyin(:,:) * gcost(:,:) - pxin(:,:) * gsint(:,:)
	122	CASE ('U') ; prot(:,:) = pyin(:,:) * gcosu(:,:) - pxin(:,:) * gsinu(:,:)
	123	CASE ('V') ; prot(:,:) = pyin(:,:) * gcosv(:,:) - pxin(:,:) * gsinv(:,:)
	124	CASE ('F') ; prot(:,:) = pyin(:,:) * gcosf(:,:) - pxin(:,:) * gsinf(:,:)
[672]	125	CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
	126	END SELECT
	127	CASE ('ij->e') ! 'ij->e' model i-j componantes to est componante
	128	SELECT CASE (cd_type)
[685]	129	CASE ('T') ; prot(:,:) = pxin(:,:) * gcost(:,:) - pyin(:,:) * gsint(:,:)
	130	CASE ('U') ; prot(:,:) = pxin(:,:) * gcosu(:,:) - pyin(:,:) * gsinu(:,:)
	131	CASE ('V') ; prot(:,:) = pxin(:,:) * gcosv(:,:) - pyin(:,:) * gsinv(:,:)
	132	CASE ('F') ; prot(:,:) = pxin(:,:) * gcosf(:,:) - pyin(:,:) * gsinf(:,:)
[672]	133	CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
	134	END SELECT
	135	CASE ('ij->n') ! 'ij->n' model i-j componantes to est componante
	136	SELECT CASE (cd_type)
[685]	137	CASE ('T') ; prot(:,:) = pyin(:,:) * gcost(:,:) + pxin(:,:) * gsint(:,:)
	138	CASE ('U') ; prot(:,:) = pyin(:,:) * gcosu(:,:) + pxin(:,:) * gsinu(:,:)
	139	CASE ('V') ; prot(:,:) = pyin(:,:) * gcosv(:,:) + pxin(:,:) * gsinv(:,:)
	140	CASE ('F') ; prot(:,:) = pyin(:,:) * gcosf(:,:) + pxin(:,:) * gsinf(:,:)
[672]	141	CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
	142	END SELECT
	143	CASE DEFAULT ; CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' )
	144	END SELECT
[3]	145
[685]	146	END SUBROUTINE rot_rep
[3]	147
	148
	149	SUBROUTINE angle
	150	!!----------------------------------------------------------------------
	151	!! * ROUTINE angle *
	152	!!
[672]	153	!! ** Purpose : Compute angles between model grid lines and the North direction
[3]	154	!!
	155	!! ** Method :
	156	!!
[672]	157	!! ** Action : Compute (gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf) arrays:
	158	!! sinus and cosinus of the angle between the north-south axe and the
	159	!! j-direction at t, u, v and f-points
[3]	160	!!
	161	!! History :
[672]	162	!! 7.0 ! 96-07 (O. Marti ) Original code
	163	!! 8.0 ! 98-06 (G. Madec )
	164	!! 8.5 ! 98-06 (G. Madec ) Free form, F90 + opt.
	165	!! 9.2 ! 07-04 (S. Masson) Add T, F points and bugfix in cos lateral boundary
[3]	166	!!----------------------------------------------------------------------
	167	INTEGER :: ji, jj ! dummy loop indices
[1884]	168	!!
[3]	169	REAL(wp) :: &
[672]	170	zlam, zphi, & ! temporary scalars
	171	zlan, zphh, & ! " "
	172	zxnpt, zynpt, znnpt, & ! x,y components and norm of the vector: T point to North Pole
	173	zxnpu, zynpu, znnpu, & ! x,y components and norm of the vector: U point to North Pole
	174	zxnpv, zynpv, znnpv, & ! x,y components and norm of the vector: V point to North Pole
	175	zxnpf, zynpf, znnpf, & ! x,y components and norm of the vector: F point to North Pole
	176	zxvvt, zyvvt, znvvt, & ! x,y components and norm of the vector: between V points below and above a T point
	177	zxffu, zyffu, znffu, & ! x,y components and norm of the vector: between F points below and above a U point
	178	zxffv, zyffv, znffv, & ! x,y components and norm of the vector: between F points left and right a V point
	179	zxuuf, zyuuf, znuuf ! x,y components and norm of the vector: between U points below and above a F point
[3]	180	!!----------------------------------------------------------------------
	181
	182	! ============================= !
	183	! Compute the cosinus and sinus !
	184	! ============================= !
	185	! (computation done on the north stereographic polar plane)
	186
[672]	187	DO jj = 2, jpjm1
[3]	188	!CDIR NOVERRCHK
	189	DO ji = fs_2, jpi ! vector opt.
	190
[672]	191	! north pole direction & modulous (at t-point)
	192	zlam = glamt(ji,jj)
	193	zphi = gphit(ji,jj)
	194	zxnpt = 0. - 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	195	zynpt = 0. - 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	196	znnpt = zxnptzxnpt + zynptzynpt
	197
[3]	198	! north pole direction & modulous (at u-point)
	199	zlam = glamu(ji,jj)
	200	zphi = gphiu(ji,jj)
	201	zxnpu = 0. - 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	202	zynpu = 0. - 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	203	znnpu = zxnpuzxnpu + zynpuzynpu
	204
	205	! north pole direction & modulous (at v-point)
	206	zlam = glamv(ji,jj)
	207	zphi = gphiv(ji,jj)
	208	zxnpv = 0. - 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	209	zynpv = 0. - 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	210	znnpv = zxnpvzxnpv + zynpvzynpv
	211
[672]	212	! north pole direction & modulous (at f-point)
	213	zlam = glamf(ji,jj)
	214	zphi = gphif(ji,jj)
	215	zxnpf = 0. - 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	216	zynpf = 0. - 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. )
	217	znnpf = zxnpfzxnpf + zynpfzynpf
	218
	219	! j-direction: v-point segment direction (around t-point)
	220	zlam = glamv(ji,jj )
	221	zphi = gphiv(ji,jj )
	222	zlan = glamv(ji,jj-1)
	223	zphh = gphiv(ji,jj-1)
	224	zxvvt = 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	225	& - 2. * COS( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	226	zyvvt = 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	227	& - 2. * SIN( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	228	znvvt = SQRT( znnpt * ( zxvvtzxvvt + zyvvtzyvvt ) )
	229	znvvt = MAX( znvvt, 1.e-14 )
	230
	231	! j-direction: f-point segment direction (around u-point)
[3]	232	zlam = glamf(ji,jj )
	233	zphi = gphif(ji,jj )
	234	zlan = glamf(ji,jj-1)
	235	zphh = gphif(ji,jj-1)
	236	zxffu = 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	237	& - 2. * COS( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	238	zyffu = 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	239	& - 2. * SIN( radzlan ) TAN( rpi/4. - rad*zphh/2. )
[672]	240	znffu = SQRT( znnpu * ( zxffuzxffu + zyffuzyffu ) )
	241	znffu = MAX( znffu, 1.e-14 )
[3]	242
[672]	243	! i-direction: f-point segment direction (around v-point)
[3]	244	zlam = glamf(ji ,jj)
	245	zphi = gphif(ji ,jj)
	246	zlan = glamf(ji-1,jj)
	247	zphh = gphif(ji-1,jj)
	248	zxffv = 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	249	& - 2. * COS( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	250	zyffv = 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	251	& - 2. * SIN( radzlan ) TAN( rpi/4. - rad*zphh/2. )
[672]	252	znffv = SQRT( znnpv * ( zxffvzxffv + zyffvzyffv ) )
	253	znffv = MAX( znffv, 1.e-14 )
[3]	254
[672]	255	! j-direction: u-point segment direction (around f-point)
	256	zlam = glamu(ji,jj+1)
	257	zphi = gphiu(ji,jj+1)
	258	zlan = glamu(ji,jj )
	259	zphh = gphiu(ji,jj )
	260	zxuuf = 2. * COS( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	261	& - 2. * COS( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	262	zyuuf = 2. * SIN( radzlam ) TAN( rpi/4. - rad*zphi/2. ) &
	263	& - 2. * SIN( radzlan ) TAN( rpi/4. - rad*zphh/2. )
	264	znuuf = SQRT( znnpf * ( zxuufzxuuf + zyuufzyuuf ) )
	265	znuuf = MAX( znuuf, 1.e-14 )
	266
[3]	267	! cosinus and sinus using scalar and vectorial products
[672]	268	gsint(ji,jj) = ( zxnptzyvvt - zynptzxvvt ) / znvvt
	269	gcost(ji,jj) = ( zxnptzxvvt + zynptzyvvt ) / znvvt
[3]	270
[672]	271	gsinu(ji,jj) = ( zxnpuzyffu - zynpuzxffu ) / znffu
	272	gcosu(ji,jj) = ( zxnpuzxffu + zynpuzyffu ) / znffu
	273
	274	gsinf(ji,jj) = ( zxnpfzyuuf - zynpfzxuuf ) / znuuf
	275	gcosf(ji,jj) = ( zxnpfzxuuf + zynpfzyuuf ) / znuuf
	276
[3]	277	! (caution, rotation of 90 degres)
[672]	278	gsinv(ji,jj) = ( zxnpvzxffv + zynpvzyffv ) / znffv
	279	gcosv(ji,jj) =-( zxnpvzyffv - zynpvzxffv ) / znffv
[3]	280
	281	END DO
	282	END DO
	283
	284	! =============== !
	285	! Geographic mesh !
	286	! =============== !
	287
[672]	288	DO jj = 2, jpjm1
[3]	289	DO ji = fs_2, jpi ! vector opt.
[672]	290	IF( MOD( ABS( glamv(ji,jj) - glamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
	291	gsint(ji,jj) = 0.
	292	gcost(ji,jj) = 1.
	293	ENDIF
	294	IF( MOD( ABS( glamf(ji,jj) - glamf(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
[3]	295	gsinu(ji,jj) = 0.
	296	gcosu(ji,jj) = 1.
	297	ENDIF
[672]	298	IF( ABS( gphif(ji,jj) - gphif(ji-1,jj) ) < 1.e-8 ) THEN
[3]	299	gsinv(ji,jj) = 0.
	300	gcosv(ji,jj) = 1.
	301	ENDIF
[672]	302	IF( MOD( ABS( glamu(ji,jj) - glamu(ji,jj+1) ), 360. ) < 1.e-8 ) THEN
	303	gsinf(ji,jj) = 0.
	304	gcosf(ji,jj) = 1.
	305	ENDIF
[3]	306	END DO
	307	END DO
	308
	309	! =========================== !
	310	! Lateral boundary conditions !
	311	! =========================== !
	312
[672]	313	! lateral boundary cond.: T-, U-, V-, F-pts, sgn
[1884]	314	CALL lbc_lnk( gcost, 'T', 1. ) ; CALL lbc_lnk( gsint, 'T', -1. )
	315	CALL lbc_lnk( gcosu, 'U', 1. ) ; CALL lbc_lnk( gsinu, 'U', -1. )
	316	CALL lbc_lnk( gcosv, 'V', 1. ) ; CALL lbc_lnk( gsinv, 'V', -1. )
	317	CALL lbc_lnk( gcosf, 'F', 1. ) ; CALL lbc_lnk( gsinf, 'F', -1. )
[3]	318
	319	END SUBROUTINE angle
	320
	321
[1884]	322	SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid, &
	323	pte, ptn )
[3]	324	!!----------------------------------------------------------------------
	325	!! * ROUTINE geo2oce *
	326	!!
	327	!! ** Purpose :
	328	!!
	329	!! ** Method : Change wind stress from geocentric to east/north
	330	!!
	331	!! History :
	332	!! ! (O. Marti) Original code
	333	!! ! 91-03 (G. Madec)
	334	!! ! 92-07 (M. Imbard)
	335	!! ! 99-11 (M. Imbard) NetCDF format with IOIPSL
	336	!! ! 00-08 (D. Ludicone) Reduced section at Bab el Mandeb
	337	!! 8.5 ! 02-06 (G. Madec) F90: Free form
[1884]	338	!! 3.0 ! 07-08 (G. Madec) geo2oce suppress lon/lat agruments
[3]	339	!!----------------------------------------------------------------------
[1884]	340	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pxx, pyy, pzz
	341	CHARACTER(len=1) , INTENT(in ) :: cgrid
	342	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pte, ptn
	343	!!
[3]	344	REAL(wp), PARAMETER :: rpi = 3.141592653E0
[88]	345	REAL(wp), PARAMETER :: rad = rpi / 180.e0
[3]	346	INTEGER :: ig !
	347	!! * Local save
[1884]	348	REAL(wp), SAVE, DIMENSION(jpi,jpj,4) :: zsinlon, zcoslon, zsinlat, zcoslat
	349	LOGICAL , SAVE, DIMENSION(4) :: linit = .FALSE.
[3]	350	!!----------------------------------------------------------------------
	351
	352	SELECT CASE( cgrid)
[1884]	353	CASE ( 'T' )
	354	ig = 1
	355	IF( .NOT. linit(ig) ) THEN
	356	zsinlon(:,:,ig) = SIN( rad * glamt(:,:) )
	357	zcoslon(:,:,ig) = COS( rad * glamt(:,:) )
	358	zsinlat(:,:,ig) = SIN( rad * gphit(:,:) )
	359	zcoslat(:,:,ig) = COS( rad * gphit(:,:) )
	360	linit(ig) = .TRUE.
	361	ENDIF
	362	CASE ( 'U' )
	363	ig = 2
	364	IF( .NOT. linit(ig) ) THEN
	365	zsinlon(:,:,ig) = SIN( rad * glamu(:,:) )
	366	zcoslon(:,:,ig) = COS( rad * glamu(:,:) )
	367	zsinlat(:,:,ig) = SIN( rad * gphiu(:,:) )
	368	zcoslat(:,:,ig) = COS( rad * gphiu(:,:) )
	369	linit(ig) = .TRUE.
	370	ENDIF
	371	CASE ( 'V' )
	372	ig = 3
	373	IF( .NOT. linit(ig) ) THEN
	374	zsinlon(:,:,ig) = SIN( rad * glamv(:,:) )
	375	zcoslon(:,:,ig) = COS( rad * glamv(:,:) )
	376	zsinlat(:,:,ig) = SIN( rad * gphiv(:,:) )
	377	zcoslat(:,:,ig) = COS( rad * gphiv(:,:) )
	378	linit(ig) = .TRUE.
	379	ENDIF
	380	CASE ( 'F' )
	381	ig = 4
	382	IF( .NOT. linit(ig) ) THEN
	383	zsinlon(:,:,ig) = SIN( rad * glamf(:,:) )
	384	zcoslon(:,:,ig) = COS( rad * glamf(:,:) )
	385	zsinlat(:,:,ig) = SIN( rad * gphif(:,:) )
	386	zcoslat(:,:,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 )
	392	END SELECT
	393
	394	pte = - zsinlon(:,:,ig) * pxx + zcoslon(:,:,ig) * pyy
	395	ptn = - zcoslon(:,:,ig) * zsinlat(:,:,ig) * pxx &
	396	- zsinlon(:,:,ig) * zsinlat(:,:,ig) * pyy &
	397	+ zcoslat(:,:,ig) * pzz
	398	!!$ ptv = zcoslon(:,:,ig) * zcoslat(:,:,ig) * pxx &
	399	!!$ + zsinlon(:,:,ig) * zcoslat(:,:,ig) * pyy &
	400	!!$ + zsinlat(:,:,ig) * pzz
	401	!
	402	END SUBROUTINE geo2oce
[3]	403
[1884]	404	SUBROUTINE oce2geo ( pte, ptn, cgrid, &
	405	pxx , pyy , pzz )
	406	!!----------------------------------------------------------------------
	407	!! * ROUTINE oce2geo *
	408	!!
	409	!! ** Purpose :
	410	!!
	411	!! ** Method : Change vector from east/north to geocentric
	412	!!
	413	!! History :
	414	!! ! (A. Caubel) oce2geo - Original code
	415	!!----------------------------------------------------------------------
	416	REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) :: pte, ptn
	417	CHARACTER(len=1) , INTENT( IN ) :: cgrid
	418	REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ) :: pxx , pyy , pzz
	419	!!
	420	REAL(wp), PARAMETER :: rpi = 3.141592653E0
	421	REAL(wp), PARAMETER :: rad = rpi / 180.e0
	422	INTEGER :: ig !
	423	!! * Local save
	424	REAL(wp), SAVE, DIMENSION(jpi,jpj,4) :: zsinlon, zcoslon, zsinlat, zcoslat
	425	LOGICAL , SAVE, DIMENSION(4) :: linit = .FALSE.
	426	!!----------------------------------------------------------------------
[3]	427
[1884]	428	SELECT CASE( cgrid)
	429	CASE ( 'T' )
	430	ig = 1
	431	IF( .NOT. linit(ig) ) THEN
	432	zsinlon(:,:,ig) = SIN( rad * glamt(:,:) )
	433	zcoslon(:,:,ig) = COS( rad * glamt(:,:) )
	434	zsinlat(:,:,ig) = SIN( rad * gphit(:,:) )
	435	zcoslat(:,:,ig) = COS( rad * gphit(:,:) )
	436	linit(ig) = .TRUE.
	437	ENDIF
	438	CASE ( 'U' )
	439	ig = 2
	440	IF( .NOT. linit(ig) ) THEN
	441	zsinlon(:,:,ig) = SIN( rad * glamu(:,:) )
	442	zcoslon(:,:,ig) = COS( rad * glamu(:,:) )
	443	zsinlat(:,:,ig) = SIN( rad * gphiu(:,:) )
	444	zcoslat(:,:,ig) = COS( rad * gphiu(:,:) )
	445	linit(ig) = .TRUE.
	446	ENDIF
	447	CASE ( 'V' )
	448	ig = 3
	449	IF( .NOT. linit(ig) ) THEN
	450	zsinlon(:,:,ig) = SIN( rad * glamv(:,:) )
	451	zcoslon(:,:,ig) = COS( rad * glamv(:,:) )
	452	zsinlat(:,:,ig) = SIN( rad * gphiv(:,:) )
	453	zcoslat(:,:,ig) = COS( rad * gphiv(:,:) )
	454	linit(ig) = .TRUE.
	455	ENDIF
	456	CASE ( 'F' )
	457	ig = 4
	458	IF( .NOT. linit(ig) ) THEN
	459	zsinlon(:,:,ig) = SIN( rad * glamf(:,:) )
	460	zcoslon(:,:,ig) = COS( rad * glamf(:,:) )
	461	zsinlat(:,:,ig) = SIN( rad * gphif(:,:) )
	462	zcoslat(:,:,ig) = COS( rad * gphif(:,:) )
	463	linit(ig) = .TRUE.
	464	ENDIF
	465	CASE default
[474]	466	WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid
	467	CALL ctl_stop( ctmp1 )
[1884]	468	END SELECT
	469
	470	pxx = - zsinlon(:,:,ig) * pte - zcoslon(:,:,ig) * zsinlat(:,:,ig) * ptn
	471	pyy = zcoslon(:,:,ig) * pte - zsinlon(:,:,ig) * zsinlat(:,:,ig) * ptn
	472	pzz = zcoslat(:,:,ig) * ptn
	473
[3]	474
[1884]	475	END SUBROUTINE oce2geo
[3]	476
	477
[672]	478	SUBROUTINE repere ( px1, py1, px2, py2, kchoix, cd_type )
[3]	479	!!----------------------------------------------------------------------
	480	!! * ROUTINE repere *
	481	!!
	482	!! ** Purpose : Change vector componantes between a geopgraphic grid
	483	!! and a stretched coordinates grid.
	484	!!
[672]	485	!! ** Method :
[3]	486	!!
	487	!! ** Action :
	488	!!
	489	!! History :
	490	!! ! 89-03 (O. Marti) original code
	491	!! ! 92-02 (M. Imbard)
	492	!! ! 93-03 (M. Guyon) symetrical conditions
	493	!! ! 98-05 (B. Blanke)
	494	!! 8.5 ! 02-08 (G. Madec) F90: Free form
	495	!!----------------------------------------------------------------------
	496	!! * Arguments
[672]	497	REAL(wp), INTENT( IN ), DIMENSION(jpi,jpj) :: &
[3]	498	px1, py1 ! two horizontal components to be rotated
[672]	499	REAL(wp), INTENT( OUT ), DIMENSION(jpi,jpj) :: &
[3]	500	px2, py2 ! the two horizontal components in the model repere
[672]	501	INTEGER, INTENT( IN ) :: &
[3]	502	kchoix ! type of transformation
	503	! = 1 change from geographic to model grid.
	504	! =-1 change from model to geographic grid
[672]	505	CHARACTER(len=1), INTENT( IN ), OPTIONAL :: cd_type ! define the nature of pt2d array grid-points
	506	!
	507	CHARACTER(len=1) :: cl_type ! define the nature of pt2d array grid-points (T point by default)
[3]	508	!!----------------------------------------------------------------------
	509
[672]	510	cl_type = 'T'
	511	IF( PRESENT(cd_type) ) cl_type = cd_type
	512	!
	513	SELECT CASE (kchoix)
	514	CASE ( 1) ! change from geographic to model grid.
[685]	515	CALL rot_rep( px1, py1, cl_type, 'en->i', px2 )
	516	CALL rot_rep( px1, py1, cl_type, 'en->j', py2 )
[672]	517	CASE (-1) ! change from model to geographic grid
[685]	518	CALL rot_rep( px1, py1, cl_type, 'ij->e', px2 )
	519	CALL rot_rep( px1, py1, cl_type, 'ij->n', py2 )
[672]	520	CASE DEFAULT ; CALL ctl_stop( 'repere: Syntax Error in the definition of kchoix (1 OR -1' )
	521	END SELECT
[3]	522
	523	END SUBROUTINE repere
	524
[1884]	525	SUBROUTINE obs_rot ( psinu, pcosu, psinv, pcosv )
	526	!!----------------------------------------------------------------------
	527	!! * ROUTINE obs_rot *
	528	!!
	529	!! ** Purpose : Copy gsinu, gcosu, gsinv and gsinv
	530	!! to input data for rotations of
	531	!! current at observation points
	532	!!
	533	!! History :
	534	!! 9.2 ! 09-02 (K. Mogensen)
	535	!!----------------------------------------------------------------------
	536	REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ):: &
	537	& psinu, pcosu, psinv, pcosv! copy of data
	538
	539	!!----------------------------------------------------------------------
	540
	541	! Initialization of gsin* and gcos* at first call
	542	! -----------------------------------------------
	543
	544	IF( lmust_init ) THEN
	545	IF(lwp) WRITE(numout,*)
	546	IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched'
	547	IF(lwp) WRITE(numout,*) ' ~~~~~~~ coordinate transformation'
	548
	549	CALL angle ! initialization of the transformation
	550	lmust_init = .FALSE.
	551
	552	ENDIF
	553
	554	psinu(:,:) = gsinu(:,:)
	555	pcosu(:,:) = gcosu(:,:)
	556	psinv(:,:) = gsinv(:,:)
	557	pcosv(:,:) = gcosv(:,:)
	558
	559	END SUBROUTINE obs_rot
	560
	561
[3]	562	!!======================================================================
	563	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: branches/TAM_V3_0/NEMO/OPA_SRC/geo2ocean.F90 @ 3294

Download in other formats: