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.

traldf_lap_blp.F90 in NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/traldf_lap_blp.F90 @ 14986

Last change on this file since 14986 was 14986, checked in by sparonuz, 3 years ago
Merge trunk -r14984:HEAD
Property svn:keywords set to `Id`
File size: 16.2 KB

Rev	Line
[5883]	1	MODULE traldf_lap_blp
[3]	2	!!==============================================================================
[5883]	3	!! * MODULE traldf_lap_blp *
[5836]	4	!! Ocean tracers: lateral diffusivity trend (laplacian and bilaplacian)
[3]	5	!!==============================================================================
[14072]	6	!! History : 3.7 ! 2014-01 (G. Madec, S. Masson) Original code, re-entrant laplacian
[2528]	7	!!----------------------------------------------------------------------
[3]	8
	9	!!----------------------------------------------------------------------
[5883]	10	!! tra_ldf_lap : tracer trend update with iso-level laplacian diffusive operator
	11	!! tra_ldf_blp : tracer trend update with iso-level or iso-neutral bilaplacian operator
[3]	12	!!----------------------------------------------------------------------
[5883]	13	USE oce ! ocean dynamics and active tracers
	14	USE dom_oce ! ocean space and time domain
[13982]	15	USE domutl, ONLY : is_tile
[5883]	16	USE ldftra ! lateral physics: eddy diffusivity
	17	USE traldf_iso ! iso-neutral lateral diffusion (standard operator) (tra_ldf_iso routine)
	18	USE traldf_triad ! iso-neutral lateral diffusion (triad operator) (tra_ldf_triad routine)
	19	USE diaptr ! poleward transport diagnostics
[7646]	20	USE diaar5 ! AR5 diagnostics
[5883]	21	USE trc_oce ! share passive tracers/Ocean variables
	22	USE zpshde ! partial step: hor. derivative (zps_hde routine)
[5836]	23	!
[5883]	24	USE in_out_manager ! I/O manager
[9019]	25	USE iom ! I/O library
[5883]	26	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	27	USE lib_mpp ! distribued memory computing library
	28	USE timing ! Timing
[3]	29
	30	IMPLICIT NONE
	31	PRIVATE
	32
[5883]	33	PUBLIC tra_ldf_lap ! called by traldf.F90
	34	PUBLIC tra_ldf_blp ! called by traldf.F90
[3]	35
[7646]	36	LOGICAL :: l_ptr ! flag to compute poleward transport
	37	LOGICAL :: l_hst ! flag to compute heat transport
	38
[3]	39	!! * Substitutions
[12377]	40	# include "do_loop_substitute.h90"
[13237]	41	# include "domzgr_substitute.h90"
[14219]	42	# include "single_precision_substitute.h90"
[3]	43	!!----------------------------------------------------------------------
[9598]	44	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[2528]	45	!! $Id$
[10068]	46	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	47	!!----------------------------------------------------------------------
	48	CONTAINS
	49
[13982]	50	SUBROUTINE tra_ldf_lap( kt, Kmm, kit000, cdtype, pahu, pahv, &
	51	& pgu , pgv , pgui, pgvi, &
	52	& pt, pt_rhs, kjpt, kpass )
	53	!!
	54	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	55	INTEGER , INTENT(in ) :: kit000 ! first time step index
	56	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	57	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	58	INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage
	59	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
	60	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s]
	61	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels
	62	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels
	63	REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pt ! before tracer fields
[14219]	64	REAL(dp), DIMENSION(:,:,:,:), INTENT(inout) :: pt_rhs ! tracer trend
[13982]	65	!!
	66	CALL tra_ldf_lap_t( kt, Kmm, kit000, cdtype, pahu, pahv, is_tile(pahu), &
	67	& pgu , pgv , is_tile(pgu) , pgui, pgvi, is_tile(pgui), &
	68	& pt, is_tile(pt), pt_rhs, is_tile(pt_rhs), kjpt, kpass )
	69	END SUBROUTINE tra_ldf_lap
	70
	71
	72	SUBROUTINE tra_ldf_lap_t( kt, Kmm, kit000, cdtype, pahu, pahv, ktah, &
	73	& pgu , pgv , ktg , pgui, pgvi, ktgi, &
	74	& pt, ktt, pt_rhs, ktt_rhs, kjpt, kpass )
[3]	75	!!----------------------------------------------------------------------
	76	!! * ROUTINE tra_ldf_lap *
[14072]	77	!!
	78	!! ** Purpose : Compute the before horizontal tracer (t & s) diffusive
[3]	79	!! trend and add it to the general trend of tracer equation.
	80	!!
	81	!! ** Method : Second order diffusive operator evaluated using before
[14072]	82	!! fields (forward time scheme). The horizontal diffusive trends of
[2528]	83	!! the tracer is given by:
[5836]	84	!! difft = 1/(e1e2te3t) { di-1[ pahu e2ue3u/e1u di(tb) ]
	85	!! + dj-1[ pahv e1v*e3v/e2v dj(tb) ] }
[12377]	86	!! Add this trend to the general tracer trend pt_rhs :
	87	!! pt_rhs = pt_rhs + difft
[3]	88	!!
[14072]	89	!! ** Action : - Update pt_rhs arrays with the before iso-level
[3]	90	!! harmonic mixing trend.
[2528]	91	!!----------------------------------------------------------------------
	92	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[5836]	93	INTEGER , INTENT(in ) :: kit000 ! first time step index
[2528]	94	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	95	INTEGER , INTENT(in ) :: kjpt ! number of tracers
[5836]	96	INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage
[12377]	97	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[13982]	98	INTEGER , INTENT(in ) :: ktah, ktg, ktgi, ktt, ktt_rhs
	99	REAL(wp), DIMENSION(A2D_T(ktah), JPK) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s]
	100	REAL(wp), DIMENSION(A2D_T(ktg), KJPT), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels
	101	REAL(wp), DIMENSION(A2D_T(ktgi), KJPT), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels
	102	REAL(wp), DIMENSION(A2D_T(ktt), JPK,KJPT), INTENT(in ) :: pt ! before tracer fields
[14219]	103	REAL(dp), DIMENSION(A2D_T(ktt_rhs),JPK,KJPT), INTENT(inout) :: pt_rhs ! tracer trend
[2715]	104	!
[13982]	105	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[14986]	106	INTEGER :: iij
[13982]	107	REAL(wp) :: zsign ! local scalars
	108	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: ztu, ztv, zaheeu, zaheev
[3]	109	!!----------------------------------------------------------------------
[3294]	110	!
[14986]	111	IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
[13982]	112	IF( kt == nit000 .AND. lwp ) THEN
	113	WRITE(numout,*)
	114	WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype, ', pass=', kpass
	115	WRITE(numout,*) '~~~~~~~~~~~ '
	116	ENDIF
	117	!
	118	l_hst = .FALSE.
	119	l_ptr = .FALSE.
	120	IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) ) l_ptr = .TRUE.
	121	IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
	122	& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
[3]	123	ENDIF
[5836]	124	!
[14986]	125	! Define pt_rhs halo points for multi-point haloes in bilaplacian case
	126	IF( nldf_tra == np_blp .AND. kpass == 1 ) THEN ; iij = nn_hls
	127	ELSE ; iij = 1
	128	ENDIF
	129
[5836]	130	! !== Initialization of metric arrays used for all tracers ==!
	131	IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0)
	132	ELSE ; zsign = -1._wp
	133	ENDIF
[13982]	134
[14986]	135	DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 ) !== First derivative (gradient) ==!
[12377]	136	zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) !!gm * umask(ji,jj,jk) pah masked!
	137	zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) !!gm * vmask(ji,jj,jk)
	138	END_3D
[5836]	139	!
	140	! ! =========== !
	141	DO jn = 1, kjpt ! tracer loop !
[14072]	142	! ! =========== !
	143	!
[14986]	144	DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 ) !== First derivative (gradient) ==!
[12377]	145	ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) )
	146	ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( pt(ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) )
	147	END_3D
[13497]	148	IF( ln_zps ) THEN ! set gradient at bottom/top ocean level
[14986]	149	DO_2D( iij, iij-1, iij, iij-1 ) ! bottom
[12377]	150	ztu(ji,jj,mbku(ji,jj)) = zaheeu(ji,jj,mbku(ji,jj)) * pgu(ji,jj,jn)
	151	ztv(ji,jj,mbkv(ji,jj)) = zaheev(ji,jj,mbkv(ji,jj)) * pgv(ji,jj,jn)
	152	END_2D
[13497]	153	IF( ln_isfcav ) THEN ! top in ocean cavities only
[14986]	154	DO_2D( iij, iij-1, iij, iij-1 )
[14072]	155	IF( miku(ji,jj) > 1 ) ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn)
	156	IF( mikv(ji,jj) > 1 ) ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn)
[12377]	157	END_2D
[2528]	158	ENDIF
[5836]	159	ENDIF
	160	!
[14986]	161	DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 ) !== Second derivative (divergence) added to the general tracer trends ==!
	162	! round brackets added to fix the order of floating point operations
	163	! needed to ensure halo 1 - halo 2 compatibility
	164	pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + ( ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) &
	165	& ) & ! bracket for halo 1 - halo 2 compatibility
	166	& + ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) &
	167	& ) & ! bracket for halo 1 - halo 2 compatibility
	168	& ) / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
[12377]	169	END_3D
[2528]	170	!
[5836]	171	! !== "Poleward" diffusive heat or salt transports ==!
	172	IF( ( kpass == 1 .AND. .NOT.ln_traldf_blp ) .OR. & !== first pass only ( laplacian) ==!
	173	( kpass == 2 .AND. ln_traldf_blp ) ) THEN !== 2nd pass only (bilaplacian) ==!
[7646]	174
	175	IF( l_ptr ) CALL dia_ptr_hst( jn, 'ldf', -ztv(:,:,:) )
	176	IF( l_hst ) CALL dia_ar5_hst( jn, 'ldf', -ztu(:,:,:), -ztv(:,:,:) )
[457]	177	ENDIF
[5836]	178	! ! ==================
	179	END DO ! end of tracer loop
	180	! ! ==================
[3294]	181	!
[13982]	182	END SUBROUTINE tra_ldf_lap_t
[5883]	183
[14072]	184
[12377]	185	SUBROUTINE tra_ldf_blp( kt, Kmm, kit000, cdtype, pahu, pahv , &
	186	& pgu , pgv , pgui, pgvi, &
	187	& pt , pt_rhs, kjpt, kldf )
[5883]	188	!!----------------------------------------------------------------------
	189	!! * ROUTINE tra_ldf_blp *
[14072]	190	!!
	191	!! ** Purpose : Compute the before lateral tracer diffusive
[5883]	192	!! trend and add it to the general trend of tracer equation.
	193	!!
	194	!! ** Method : The lateral diffusive trends is provided by a bilaplacian
	195	!! operator applied to before field (forward in time).
	196	!! It is computed by two successive calls to laplacian routine
	197	!!
	198	!! ** Action : pta updated with the before rotated bilaplacian diffusion
	199	!!----------------------------------------------------------------------
	200	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	201	INTEGER , INTENT(in ) :: kit000 ! first time step index
	202	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	203	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	204	INTEGER , INTENT(in ) :: kldf ! type of operator used
[12377]	205	INTEGER , INTENT(in ) :: Kmm ! ocean time level indices
[5883]	206	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s]
	207	REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels
	208	REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels
[12377]	209	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! before and now tracer fields
[14219]	210	REAL(dp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend
[5883]	211	!
	212	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[14219]	213	REAL(dp), DIMENSION(A2D(nn_hls),jpk,kjpt) :: zlap ! laplacian at t-point
[13982]	214	REAL(wp), DIMENSION(A2D(nn_hls), kjpt) :: zglu, zglv ! bottom GRADh of the laplacian (u- and v-points)
	215	REAL(wp), DIMENSION(A2D(nn_hls), kjpt) :: zgui, zgvi ! top GRADh of the laplacian (u- and v-points)
[5883]	216	!!---------------------------------------------------------------------
	217	!
[14986]	218	IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
[13982]	219	IF( kt == kit000 .AND. lwp ) THEN
	220	WRITE(numout,*)
	221	SELECT CASE ( kldf )
	222	CASE ( np_blp ) ; WRITE(numout,*) 'tra_ldf_blp : iso-level bilaplacian operator on ', cdtype
	223	CASE ( np_blp_i ) ; WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (Standard)'
	224	CASE ( np_blp_it ) ; WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (triad)'
	225	END SELECT
	226	WRITE(numout,*) '~~~~~~~~~~~'
	227	ENDIF
[5883]	228	ENDIF
	229
	230	zlap(:,:,:,:) = 0._wp
	231	!
[12377]	232	SELECT CASE ( kldf ) !== 1st laplacian applied to pt (output in zlap) ==!
[5883]	233	!
	234	CASE ( np_blp ) ! iso-level bilaplacian
[12377]	235	CALL tra_ldf_lap ( kt, Kmm, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, zlap, kjpt, 1 )
[5883]	236	CASE ( np_blp_i ) ! rotated bilaplacian : standard operator (Madec)
[12377]	237	CALL tra_ldf_iso ( kt, Kmm, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 )
[5883]	238	CASE ( np_blp_it ) ! rotated bilaplacian : triad operator (griffies)
[12377]	239	CALL tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 )
[5883]	240	END SELECT
	241	!
[14986]	242	IF (nn_hls==1) CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1.0_wp ) ! Lateral boundary conditions (unchanged sign)
[14072]	243	! ! Partial top/bottom cell: GRADh( zlap )
[14286]	244	IF( ln_isfcav .AND. ln_zps ) THEN ; CALL zps_hde_isf( kt, Kmm, kjpt, zlap, zglu, zglv, zgui, zgvi ) ! both top & bottom
	245	ELSEIF( ln_zps ) THEN ; CALL zps_hde ( kt, Kmm, kjpt, zlap, zglu, zglv ) ! only bottom
[5883]	246	ENDIF
	247	!
[12377]	248	SELECT CASE ( kldf ) !== 2nd laplacian applied to zlap (output in pt_rhs) ==!
[5883]	249	!
	250	CASE ( np_blp ) ! iso-level bilaplacian
[14219]	251	CALL tra_ldf_lap ( kt, Kmm, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, CASTWP(zlap), pt_rhs, kjpt, 2 )
[5883]	252	CASE ( np_blp_i ) ! rotated bilaplacian : standard operator (Madec)
[14219]	253	CALL tra_ldf_iso ( kt, Kmm, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, CASTWP(zlap), pt , pt_rhs, kjpt, 2 )
[5883]	254	CASE ( np_blp_it ) ! rotated bilaplacian : triad operator (griffies)
[14219]	255	CALL tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, CASTWP(zlap), pt , pt_rhs, kjpt, 2 )
[5883]	256	END SELECT
	257	!
	258	END SUBROUTINE tra_ldf_blp
	259
[3]	260	!!==============================================================================
[5883]	261	END MODULE traldf_lap_blp

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

Download in other formats: