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.

zpshde.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 10 years ago
Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.
Property svn:keywords set to `Id`
File size: 12.0 KB

Rev	Line
[3]	1	MODULE zpshde
[2528]	2	!!======================================================================
[3]	3	!! * MODULE zpshde *
[2528]	4	!! z-coordinate + partial step : Horizontal Derivative at ocean bottom level
	5	!!======================================================================
	6	!! History : OPA ! 2002-04 (A. Bozec) Original code
	7	!! NEMO 1.0 ! 2002-08 (G. Madec E. Durand) Optimization and Free form
	8	!! - ! 2004-03 (C. Ethe) adapted for passive tracers
	9	!! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA
	10	!!======================================================================
[457]	11
[3]	12	!!----------------------------------------------------------------------
	13	!! zps_hde : Horizontal DErivative of T, S and rd at the last
	14	!! ocean level (Z-coord. with Partial Steps)
	15	!!----------------------------------------------------------------------
[2528]	16	USE oce ! ocean: dynamics and tracers variables
	17	USE dom_oce ! domain: ocean variables
[3]	18	USE phycst ! physical constants
[2528]	19	USE eosbn2 ! ocean equation of state
[3]	20	USE in_out_manager ! I/O manager
	21	USE lbclnk ! lateral boundary conditions (or mpp link)
[2715]	22	USE lib_mpp ! MPP library
[3]	23
	24	IMPLICIT NONE
	25	PRIVATE
	26
[2528]	27	PUBLIC zps_hde ! routine called by step.F90
[3]	28
[3211]	29	!! * Control permutation of array indices
	30	# include "oce_ftrans.h90"
	31	# include "dom_oce_ftrans.h90"
	32
[3]	33	!! * Substitutions
	34	# include "domzgr_substitute.h90"
	35	# include "vectopt_loop_substitute.h90"
	36	!!----------------------------------------------------------------------
[2528]	37	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
	38	!! $Id$
	39	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[247]	40	!!----------------------------------------------------------------------
[3]	41	CONTAINS
	42
[2528]	43	SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv, &
	44	prd, pgru, pgrv )
[3]	45	!!----------------------------------------------------------------------
	46	!! * ROUTINE zps_hde *
	47	!!
[2528]	48	!! ** Purpose : Compute the horizontal derivative of T, S and rho
[3]	49	!! at u- and v-points with a linear interpolation for z-coordinate
	50	!! with partial steps.
	51	!!
	52	!! ** Method : In z-coord with partial steps, scale factors on last
	53	!! levels are different for each grid point, so that T, S and rd
	54	!! points are not at the same depth as in z-coord. To have horizontal
	55	!! gradients again, we interpolate T and S at the good depth :
	56	!! Linear interpolation of T, S
	57	!! Computation of di(tb) and dj(tb) by vertical interpolation:
	58	!! di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~
	59	!! dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~
	60	!! This formulation computes the two cases:
	61	!! CASE 1 CASE 2
	62	!! k-1 ___ ___________ k-1 ___ ___________
	63	!! Ti T~ T~ Ti+1
	64	!! _____ _____
	65	!! k \| \|Ti+1 k Ti \| \|
	66	!! \| \|____ ____\| \|
	67	!! ___ \| \| \| ___ \| \| \|
	68	!!
	69	!! case 1-> e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then
	70	!! t~ = t(i+1,j ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1)
	71	!! ( t~ = t(i ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1) )
	72	!! or
	73	!! case 2-> e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then
	74	!! t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i )
	75	!! ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) )
	76	!! Idem for di(s) and dj(s)
	77	!!
[87]	78	!! For rho, we call eos_insitu_2d which will compute rd~(t~,s~) at
[3]	79	!! the good depth zh from interpolated T and S for the different
	80	!! formulation of the equation of state (eos).
	81	!! Gradient formulation for rho :
[2528]	82	!! di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~
[3]	83	!!
[2528]	84	!! ** Action : - pgtu, pgtv: horizontal gradient of tracer at u- & v-points
	85	!! - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points
	86	!!----------------------------------------------------------------------
[2715]	87	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	88	USE wrk_nemo, ONLY: zri => wrk_2d_1 , zrj => wrk_2d_2 ! interpolated value of rd
	89	USE wrk_nemo, ONLY: zhi => wrk_2d_3 , zhj => wrk_2d_4 ! depth of interpolation for eos2d
	90	!
[2528]	91	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	92	INTEGER , INTENT(in ) :: kjpt ! number of tracers
[3211]	93
	94	!! DCSE_NEMO: This style defeats ftrans
	95	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pta ! 4D tracers fields
	96	!FTRANS pta :I :I :z :
[4409]	97	REAL(wp), INTENT(in ) :: pta(jpi,jpj,jpkorig,kjpt) ! 4D tracers fields
[3211]	98
[2528]	99	REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts
[3211]	100
	101	! REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields
	102	!FTRANS prd :I :I :z
[4409]	103	REAL(wp), INTENT(in ), OPTIONAL :: prd(jpi,jpj,jpkorig) ! 3D density anomaly fields
[3211]	104
[2528]	105	REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad. of prd at u- & v-pts
[2715]	106	!
[2528]	107	INTEGER :: ji, jj, jn ! Dummy loop indices
	108	INTEGER :: iku, ikv, ikum1, ikvm1 ! partial step level (ocean bottom level) at u- and v-points
	109	REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv ! temporary scalars
[2715]	110	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zti, ztj ! interpolated value of tracer
[3]	111	!!----------------------------------------------------------------------
	112
[2715]	113	IF( wrk_in_use(2, 1,2,3,4) ) THEN
	114	CALL ctl_stop('zps_hde: requested workspace arrays unavailable') ; RETURN
	115	END IF
	116
	117	! Allocate workspaces whose dimension is > jpk
	118	ALLOCATE( zti(jpi,jpj,kjpt) )
	119	ALLOCATE( ztj(jpi,jpj,kjpt) )
	120
[2528]	121	DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==!
	122	!
	123	# if defined key_vectopt_loop
	124	jj = 1
	125	DO ji = 1, jpij-jpi ! vector opt. (forced unrolled)
	126	# else
[3]	127	DO jj = 1, jpjm1
[2528]	128	DO ji = 1, jpim1
	129	# endif
[2569]	130	iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points
	131	ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1
[2528]	132	ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku)
	133	ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv)
	134	!
	135	! i- direction
	136	IF( ze3wu >= 0._wp ) THEN ! case 1
	137	zmaxu = ze3wu / fse3w(ji+1,jj,iku)
	138	! interpolated values of tracers
	139	zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
	140	! gradient of tracers
[3211]	141	#if defined key_z_first
	142	pgtu(ji,jj,jn) = umask_1(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
	143	#else
[2528]	144	pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
[3211]	145	#endif
[2528]	146	ELSE ! case 2
	147	zmaxu = -ze3wu / fse3w(ji,jj,iku)
	148	! interpolated values of tracers
	149	zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
	150	! gradient of tracers
[3211]	151	#if defined key_z_first
	152	pgtu(ji,jj,jn) = umask_1(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
	153	#else
[2528]	154	pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
[3211]	155	#endif
[2528]	156	ENDIF
	157	!
	158	! j- direction
	159	IF( ze3wv >= 0._wp ) THEN ! case 1
	160	zmaxv = ze3wv / fse3w(ji,jj+1,ikv)
	161	! interpolated values of tracers
	162	ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
	163	! gradient of tracers
[3211]	164	#if defined key_z_first
	165	pgtv(ji,jj,jn) = vmask_1(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
	166	#else
[2528]	167	pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
[3211]	168	#endif
[2528]	169	ELSE ! case 2
	170	zmaxv = -ze3wv / fse3w(ji,jj,ikv)
	171	! interpolated values of tracers
	172	ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
	173	! gradient of tracers
[3211]	174	#if defined key_z_first
	175	pgtv(ji,jj,jn) = vmask_1(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
	176	#else
[2528]	177	pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
[3211]	178	#endif
[2528]	179	ENDIF
	180	# if ! defined key_vectopt_loop
[3]	181	END DO
[2528]	182	# endif
[3]	183	END DO
[2528]	184	CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. ) ! Lateral boundary cond.
	185	!
	186	END DO
[3]	187
[2528]	188	! horizontal derivative of density anomalies (rd)
	189	IF( PRESENT( prd ) ) THEN ! depth of the partial step level
[789]	190	# if defined key_vectopt_loop
[3]	191	jj = 1
	192	DO ji = 1, jpij-jpi ! vector opt. (forced unrolled)
	193	# else
[2528]	194	DO jj = 1, jpjm1
	195	DO ji = 1, jpim1
[3]	196	# endif
[2528]	197	iku = mbku(ji,jj)
	198	ikv = mbkv(ji,jj)
	199	ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku)
	200	ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv)
	201	IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = fsdept(ji ,jj,iku) ! i-direction: case 1
	202	ELSE ; zhi(ji,jj) = fsdept(ji+1,jj,iku) ! - - case 2
	203	ENDIF
	204	IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = fsdept(ji,jj ,ikv) ! j-direction: case 1
	205	ELSE ; zhj(ji,jj) = fsdept(ji,jj+1,ikv) ! - - case 2
	206	ENDIF
[789]	207	# if ! defined key_vectopt_loop
[2528]	208	END DO
	209	# endif
[3]	210	END DO
	211
[2528]	212	! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
	213	! step and store it in zri, zrj for each case
	214	CALL eos( zti, zhi, zri ) ; CALL eos( ztj, zhj, zrj )
[3]	215
[2528]	216	! Gradient of density at the last level
[789]	217	# if defined key_vectopt_loop
[3]	218	jj = 1
	219	DO ji = 1, jpij-jpi ! vector opt. (forced unrolled)
	220	# else
[2528]	221	DO jj = 1, jpjm1
	222	DO ji = 1, jpim1
[3]	223	# endif
[2528]	224	iku = mbku(ji,jj)
	225	ikv = mbkv(ji,jj)
	226	ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku)
	227	ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv)
	228	IF( ze3wu >= 0._wp ) THEN ; pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji ,jj) - prd(ji,jj,iku) ) ! i: 1
	229	ELSE ; pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) ! i: 2
	230	ENDIF
	231	IF( ze3wv >= 0._wp ) THEN ; pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1
	232	ELSE ; pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) ! j: 2
	233	ENDIF
[789]	234	# if ! defined key_vectopt_loop
[2528]	235	END DO
	236	# endif
[3]	237	END DO
[2528]	238	CALL lbc_lnk( pgru , 'U', -1. ) ; CALL lbc_lnk( pgrv , 'V', -1. ) ! Lateral boundary conditions
	239	!
	240	END IF
	241	!
[2715]	242	IF( wrk_not_released(2, 1,2,3,4) ) CALL ctl_stop('zps_hde: failed to release workspace arrays')
	243	!
	244	DEALLOCATE( zti )
	245	DEALLOCATE( ztj )
	246	!
[3]	247	END SUBROUTINE zps_hde
	248
	249	!!======================================================================
	250	END MODULE zpshde

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

Download in other formats: