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

Context Navigation

source: branches/2013/dev_r3987_UKMO4_OBS/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90 @ 4113

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

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

Download in other formats: