Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

trazdf_imp.F90 @ 9616

Last change on this file since 9616 was 9176, checked in by andmirek, 6 years ago
#2001: OMP directives
File size: 12.5 KB

Rev	Line
[3]	1	MODULE trazdf_imp
[1438]	2	!!======================================================================
[457]	3	!! * MODULE trazdf_imp *
[2528]	4	!! Ocean tracers: vertical component of the tracer mixing trend
[1438]	5	!!======================================================================
	6	!! History : OPA ! 1990-10 (B. Blanke) Original code
	7	!! 7.0 ! 1991-11 (G. Madec)
	8	!! ! 1992-06 (M. Imbard) correction on tracer trend loops
	9	!! ! 1996-01 (G. Madec) statement function for e3
	10	!! ! 1997-05 (G. Madec) vertical component of isopycnal
	11	!! ! 1997-07 (G. Madec) geopotential diffusion in s-coord
	12	!! ! 2000-08 (G. Madec) double diffusive mixing
	13	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
	14	!! 2.0 ! 2006-11 (G. Madec) New step reorganisation
	15	!! 3.2 ! 2009-03 (G. Madec) heat and salt content trends
[2528]	16	!! 3.3 ! 2010-06 (C. Ethe, G. Madec) Merge TRA-TRC
[2602]	17	!! - ! 2011-02 (A. Coward, C. Ethe, G. Madec) improvment of surface boundary condition
[3]	18	!!----------------------------------------------------------------------
[1438]	19
	20	!!----------------------------------------------------------------------
[457]	21	!! tra_zdf_imp : Update the tracer trend with the diagonal vertical
	22	!! part of the mixing tensor.
[3]	23	!!----------------------------------------------------------------------
[457]	24	USE oce ! ocean dynamics and tracers variables
	25	USE dom_oce ! ocean space and time domain variables
	26	USE zdf_oce ! ocean vertical physics variables
[2528]	27	USE trc_oce ! share passive tracers/ocean variables
	28	USE domvvl ! variable volume
[216]	29	USE ldftra_oce ! ocean active tracers: lateral physics
[2528]	30	USE ldftra ! lateral mixing type
[457]	31	USE ldfslp ! lateral physics: slope of diffusion
	32	USE zdfddm ! ocean vertical physics: double diffusion
[2528]	33	USE traldf_iso_grif ! active tracers: Griffies operator
[3]	34	USE in_out_manager ! I/O manager
[457]	35	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[2715]	36	USE lib_mpp ! MPP library
[3294]	37	USE wrk_nemo ! Memory Allocation
	38	USE timing ! Timing
[3]	39
	40	IMPLICIT NONE
	41	PRIVATE
	42
[1438]	43	PUBLIC tra_zdf_imp ! routine called by step.F90
[3]	44
[2602]	45	REAL(wp) :: r_vvl ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise
	46
[3]	47	!! * Substitutions
	48	# include "domzgr_substitute.h90"
[457]	49	# include "ldftra_substitute.h90"
[3]	50	# include "zdfddm_substitute.h90"
[457]	51	# include "vectopt_loop_substitute.h90"
[3]	52	!!----------------------------------------------------------------------
[2528]	53	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1156]	54	!! $Id$
[2528]	55	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[457]	56	!!----------------------------------------------------------------------
[3]	57	CONTAINS
[2528]	58
[3294]	59	SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt )
[3]	60	!!----------------------------------------------------------------------
	61	!! * ROUTINE tra_zdf_imp *
	62	!!
[2602]	63	!! ** Purpose : Compute the after tracer through a implicit computation
	64	!! of the vertical tracer diffusion (including the vertical component
	65	!! of lateral mixing (only for 2nd order operator, for fourth order
	66	!! it is already computed and add to the general trend in traldf)
[3]	67	!!
[2602]	68	!! ** Method : The vertical diffusion of the tracer t is given by:
	69	!! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
[457]	70	!! It is computed using a backward time scheme (t=ta).
[2602]	71	!! If lk_zdfddm=T, use avs for salinity or for passive tracers
[3]	72	!! Surface and bottom boundary conditions: no diffusive flux on
	73	!! both tracers (bottom, applied through the masked field avt).
[2602]	74	!! If iso-neutral mixing, add to avt the contribution due to lateral mixing.
[3]	75	!!
[2602]	76	!! ** Action : - pta becomes the after tracer
[3]	77	!!---------------------------------------------------------------------
[2715]	78	USE oce , ONLY: zwd => ua , zws => va ! (ua,va) used as 3D workspace
	79	!
[2528]	80	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[3294]	81	INTEGER , INTENT(in ) :: kit000 ! first time step index
[2528]	82	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	83	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	84	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
	85	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields
	86	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[2715]	87	!
	88	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	89	REAL(wp) :: zrhs, ze3tb, ze3tn, ze3ta ! local scalars
[3294]	90	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwt
[3]	91	!!---------------------------------------------------------------------
[3294]	92	!
	93	IF( nn_timing == 1 ) CALL timing_start('tra_zdf_imp')
	94	!
	95	CALL wrk_alloc( jpi, jpj, jpk, zwi, zwt )
	96	!
	97	IF( kt == kit000 ) THEN
[457]	98	IF(lwp)WRITE(numout,*)
[2528]	99	IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype
[457]	100	IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
[2602]	101	!
	102	IF( lk_vvl ) THEN ; r_vvl = 1._wp ! Variable volume indicator
	103	ELSE ; r_vvl = 0._wp
	104	ENDIF
[457]	105	ENDIF
[2528]	106	!
[2602]	107	! ! ============= !
	108	DO jn = 1, kjpt ! tracer loop !
	109	! ! ============= !
[2528]	110	!
	111	! Matrix construction
[2602]	112	! --------------------
	113	! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
	114	!
[2715]	115	IF( ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. lk_zdfddm ) ) ) .OR. &
[2602]	116	& ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
	117	!
	118	! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
	119	IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN ; zwt(:,:,2:jpk) = avt (:,:,2:jpk)
	120	ELSE ; zwt(:,:,2:jpk) = fsavs(:,:,2:jpk)
[2528]	121	ENDIF
[4990]	122	DO jj=1, jpj
	123	DO ji=1, jpi
[5120]	124	zwt(ji,jj,1) = 0._wp
[4990]	125	END DO
	126	END DO
	127	!
[2528]	128	#if defined key_ldfslp
[2602]	129	! isoneutral diffusion: add the contribution
	130	IF( ln_traldf_grif ) THEN ! Griffies isoneutral diff
[9176]	131	!$OMP PARALLEL DO
[2528]	132	DO jk = 2, jpkm1
	133	DO jj = 2, jpjm1
	134	DO ji = fs_2, fs_jpim1 ! vector opt.
[2602]	135	zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
[2528]	136	END DO
	137	END DO
	138	END DO
[2602]	139	ELSE IF( l_traldf_rot ) THEN ! standard isoneutral diff
[9176]	140	!$OMP PARALLEL DO
[2528]	141	DO jk = 2, jpkm1
	142	DO jj = 2, jpjm1
	143	DO ji = fs_2, fs_jpim1 ! vector opt.
[2602]	144	zwt(ji,jj,jk) = zwt(ji,jj,jk) + fsahtw(ji,jj,jk) &
	145	& * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
	146	& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )
[2528]	147	END DO
	148	END DO
	149	END DO
	150	ENDIF
[592]	151	#endif
[2602]	152	! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked)
[9176]	153	!$OMP PARALLEL DO PRIVATE(ze3ta, ze3tn)
[2528]	154	DO jk = 1, jpkm1
	155	DO jj = 2, jpjm1
	156	DO ji = fs_2, fs_jpim1 ! vector opt.
[2602]	157	ze3ta = ( 1. - r_vvl ) + r_vvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point
	158	ze3tn = r_vvl + ( 1. - r_vvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point
[2528]	159	zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) )
	160	zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
	161	zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
	162	END DO
	163	END DO
[3]	164	END DO
[2528]	165	!
[2602]	166	!! Matrix inversion from the first level
	167	!!----------------------------------------------------------------------
	168	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	169	!
	170	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	171	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	172	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	173	! ( ... )( ... ) ( ... )
	174	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	175	!
	176	! m is decomposed in the product of an upper and lower triangular matrix.
	177	! The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
	178	! Suffices i,s and d indicate "inferior" (below diagonal), diagonal
	179	! and "superior" (above diagonal) components of the tridiagonal system.
	180	! The solution will be in the 4d array pta.
	181	! The 3d array zwt is used as a work space array.
	182	! En route to the solution pta is used a to evaluate the rhs and then
	183	! used as a work space array: its value is modified.
	184	!
[2528]	185	! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k)
[2602]	186	! done once for all passive tracers (so included in the IF instruction)
[2528]	187	DO jj = 2, jpjm1
	188	DO ji = fs_2, fs_jpim1
[5120]	189	zwt(ji,jj,1) = zwd(ji,jj,1)
	190	END DO
	191	END DO
[9176]	192
	193	!$OMP PARALLEL
[5120]	194	DO jk = 2, jpkm1
[9176]	195	!$OMP DO
[5120]	196	DO jj = 2, jpjm1
	197	DO ji = fs_2, fs_jpim1
[4990]	198	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
[2528]	199	END DO
	200	END DO
[9176]	201	!$OMP END DO
[2528]	202	END DO
[9176]	203	!$OMP END PARALLEL
[2528]	204	!
	205	END IF
[2602]	206	!
[2528]	207	! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1
[9176]	208	!$OMP PARALLEL DO PRIVATE(ze3tb, ze3tn)
[457]	209	DO jj = 2, jpjm1
	210	DO ji = fs_2, fs_jpim1
[5120]	211	ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1)
	212	ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1)
	213	pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) &
	214	& + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
	215	END DO
	216	END DO
[9176]	217
	218	!$OMP PARALLEL
[5120]	219	DO jk = 2, jpkm1
[9176]	220	!$OMP DO PRIVATE(ze3tb, ze3tn, zrhs)
[5120]	221	DO jj = 2, jpjm1
	222	DO ji = fs_2, fs_jpim1
[2602]	223	ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk)
	224	ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t (ji,jj,jk)
[2528]	225	zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn) ! zrhs=right hand side
	226	pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
	227	END DO
[3]	228	END DO
[9176]	229	!$OMP END DO
[3]	230	END DO
[9176]	231	!$OMP END PARALLEL
[2602]	232	! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer)
[9176]	233	!$OMP PARALLEL
	234	!$OMP DO
[457]	235	DO jj = 2, jpjm1
	236	DO ji = fs_2, fs_jpim1
[2528]	237	pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
[5120]	238	END DO
	239	END DO
[9176]	240
[5120]	241	DO jk = jpk-2, 1, -1
[9176]	242	!$OMP DO
[5120]	243	DO jj = 2, jpjm1
	244	DO ji = fs_2, fs_jpim1
[2528]	245	pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) &
	246	& / zwt(ji,jj,jk) * tmask(ji,jj,jk)
	247	END DO
[457]	248	END DO
[9176]	249	!$OMP END DO
[457]	250	END DO
[9176]	251	!$OMP END PARALLEL
[2602]	252	! ! ================= !
	253	END DO ! end tracer loop !
	254	! ! ================= !
[1438]	255	!
[3294]	256	CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt )
[2715]	257	!
[3294]	258	IF( nn_timing == 1 ) CALL timing_stop('tra_zdf_imp')
	259	!
[3]	260	END SUBROUTINE tra_zdf_imp
	261
	262	!!==============================================================================
	263	END MODULE trazdf_imp

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 9616

Download in other formats: