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

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90 @ 10395

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

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

Download in other formats: