Context Navigation

traldf_iso_grif.F90 @ 13320

Last change on this file since 13320 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 22.8 KB

Rev	Line
[2205]	1	MODULE traldf_iso_grif
[2371]	2	!!======================================================================
[2205]	3	!! * MODULE traldf_iso_grif *
[2371]	4	!! Ocean tracers: horizontal component of the lateral tracer mixing trend
	5	!!======================================================================
[3294]	6	!! History : 3.3 ! 2010-10 (G. Nurser, C. Harris, G. Madec)
[2371]	7	!! ! Griffies operator version adapted from traldf_iso.F90
[2205]	8	!!----------------------------------------------------------------------
	9	#if defined key_ldfslp \|\| defined key_esopa
	10	!!----------------------------------------------------------------------
	11	!! 'key_ldfslp' slope of the lateral diffusive direction
	12	!!----------------------------------------------------------------------
[3294]	13	!! tra_ldf_iso_grif : update the tracer trend with the horizontal component
	14	!! of the Griffies iso-neutral laplacian operator
[2205]	15	!!----------------------------------------------------------------------
	16	USE oce ! ocean dynamics and active tracers
	17	USE dom_oce ! ocean space and time domain
[2715]	18	USE phycst ! physical constants
[2454]	19	USE trc_oce ! share passive tracers/Ocean variables
	20	USE zdf_oce ! ocean vertical physics
[2205]	21	USE ldftra_oce ! ocean active tracers: lateral physics
	22	USE ldfslp ! iso-neutral slopes
	23	USE diaptr ! poleward transport diagnostics
[2454]	24	USE in_out_manager ! I/O manager
	25	USE iom ! I/O library
[2371]	26	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[2715]	27	USE lib_mpp ! MPP library
[3294]	28	USE wrk_nemo ! Memory Allocation
	29	USE timing ! Timing
[2205]	30
[3294]	31
[2205]	32	IMPLICIT NONE
	33	PRIVATE
	34
[2715]	35	PUBLIC tra_ldf_iso_grif ! routine called by traldf.F90
[2205]	36
[2715]	37	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: psix_eiv, psiy_eiv !: eiv stream function (diag only)
	38	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: ah_wslp2 !: aeiv*w-slope^2
[3294]	39	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zdkt3d !: vertical tracer gradient at 2 levels
[2371]	40
[2205]	41	!! * Substitutions
	42	# include "domzgr_substitute.h90"
	43	# include "ldftra_substitute.h90"
[2371]	44	# include "vectopt_loop_substitute.h90"
[2205]	45	# include "ldfeiv_substitute.h90"
	46	!!----------------------------------------------------------------------
[2287]	47	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
	48	!! $Id$
[2399]	49	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[2205]	50	!!----------------------------------------------------------------------
	51	CONTAINS
	52
[3294]	53	SUBROUTINE tra_ldf_iso_grif( kt, kit000, cdtype, pgu, pgv, &
[2399]	54	& ptb, pta, kjpt, pahtb0 )
[2450]	55	!!----------------------------------------------------------------------
	56	!! * ROUTINE tra_ldf_iso_grif *
	57	!!
[3294]	58	!! ** Purpose : Compute the before horizontal tracer (t & s) diffusive
	59	!! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
[2450]	60	!! add it to the general trend of tracer equation.
	61	!!
[3294]	62	!! ** Method : The horizontal component of the lateral diffusive trends
[2450]	63	!! is provided by a 2nd order operator rotated along neural or geopo-
	64	!! tential surfaces to which an eddy induced advection can be added
	65	!! It is computed using before fields (forward in time) and isopyc-
	66	!! nal or geopotential slopes computed in routine ldfslp.
	67	!!
	68	!! 1st part : masked horizontal derivative of T ( di[ t ] )
	69	!! ======== with partial cell update if ln_zps=T.
	70	!!
	71	!! 2nd part : horizontal fluxes of the lateral mixing operator
[3294]	72	!! ========
[2450]	73	!! zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ]
	74	!! - aht e2u*uslp dk[ mi(mk(tb)) ]
	75	!! zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ]
	76	!! - aht e2u*vslp dk[ mj(mk(tb)) ]
	77	!! take the horizontal divergence of the fluxes:
	78	!! difft = 1/(e1te2te3t) { di-1[ zftu ] + dj-1[ zftv ] }
	79	!! Add this trend to the general trend (ta,sa):
	80	!! ta = ta + difft
	81	!!
	82	!! 3rd part: vertical trends of the lateral mixing operator
	83	!! ======== (excluding the vertical flux proportional to dk[t] )
	84	!! vertical fluxes associated with the rotated lateral mixing:
	85	!! zftw =-aht { e2t*wslpi di[ mi(mk(tb)) ]
	86	!! + e1t*wslpj dj[ mj(mk(tb)) ] }
	87	!! take the horizontal divergence of the fluxes:
	88	!! difft = 1/(e1te2te3t) dk[ zftw ]
	89	!! Add this trend to the general trend (ta,sa):
	90	!! pta = pta + difft
	91	!!
	92	!! ** Action : Update pta arrays with the before rotated diffusion
	93	!!----------------------------------------------------------------------
[2715]	94	USE oce , ONLY: zftu => ua , zftv => va ! (ua,va) used as 3D workspace
	95	!
[2450]	96	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[3294]	97	INTEGER , INTENT(in ) :: kit000 ! first time step index
[2450]	98	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	99	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	100	REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels
	101	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields
[3294]	102	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[2450]	103	REAL(wp) , INTENT(in ) :: pahtb0 ! background diffusion coef
[2715]	104	!
[2450]	105	INTEGER :: ji, jj, jk,jn ! dummy loop indices
	106	INTEGER :: ip,jp,kp ! dummy loop indices
	107	INTEGER :: ierr ! temporary integer
	108	REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars
	109	REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - -
	110	REAL(wp) :: zcoef0, zbtr ! - -
[2371]	111	!
[2454]	112	REAL(wp) :: zslope_skew, zslope_iso, zslope2, zbu, zbv
	113	REAL(wp) :: ze1ur, zdxt, ze2vr, ze3wr, zdyt, zdzt
	114	REAL(wp) :: zah, zah_slp, zaei_slp
[3294]	115	REAL(wp), POINTER, DIMENSION(:,: ) :: z2d
	116	REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, ztfw
	117	REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d ! 3D workspace
[2205]	118	!!----------------------------------------------------------------------
[3294]	119	!
	120	IF( nn_timing == 1 ) CALL timing_start('tra_ldf_iso_grif')
	121	!
	122	CALL wrk_alloc( jpi, jpj, z2d )
	123	CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw )
	124	!
[2205]	125
[3294]	126	IF( kt == kit000 .AND. .NOT.ALLOCATED(ah_wslp2) ) THEN
[2450]	127	IF(lwp) WRITE(numout,*)
	128	IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator on ', cdtype
	129	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[11101]	130	IF(lwp .AND. lflush) CALL flush(numout)
[3294]	131	ALLOCATE( ah_wslp2(jpi,jpj,jpk) , zdkt3d(jpi,jpj,0:1), STAT=ierr )
[2715]	132	IF( lk_mpp ) CALL mpp_sum ( ierr )
	133	IF( ierr > 0 ) CALL ctl_stop('STOP', 'tra_ldf_iso_grif: unable to allocate arrays')
[2450]	134	IF( ln_traldf_gdia ) THEN
[3294]	135	IF (.NOT. ALLOCATED(psix_eiv))THEN
	136	ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr )
	137	IF( lk_mpp ) CALL mpp_sum ( ierr )
	138	IF( ierr > 0 ) CALL ctl_stop('STOP', 'tra_ldf_iso_grif: unable to allocate diagnostics')
	139	ENDIF
[2450]	140	ENDIF
	141	ENDIF
[2371]	142
[2205]	143	!!----------------------------------------------------------------------
[3294]	144	!! 0 - calculate ah_wslp2, psix_eiv, psiy_eiv
[2371]	145	!!----------------------------------------------------------------------
[2205]	146
[3294]	147	!!gm Future development: consider using Ah defined at T-points and attached to the 4 t-point triads
[2205]	148
[2454]	149	ah_wslp2(:,:,:) = 0._wp
[2450]	150	IF( ln_traldf_gdia ) THEN
[2454]	151	psix_eiv(:,:,:) = 0._wp
	152	psiy_eiv(:,:,:) = 0._wp
[2450]	153	ENDIF
[2205]	154
[2454]	155	DO ip = 0, 1
	156	DO kp = 0, 1
	157	DO jk = 1, jpkm1
	158	DO jj = 1, jpjm1
	159	DO ji = 1, fs_jpim1
[3294]	160	ze1ur = 1._wp / e1u(ji,jj)
[2454]	161	ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
	162	zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
[3294]	163	zah = fsahtu(ji,jj,jk) ! fsaht(ji+ip,jj,jk)
[2450]	164	zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
[3294]	165	! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
	166	! (do this by adding gradient of depth)
	167	zslope2 = zslope_skew + ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * ze1ur * umask(ji,jj,jk+kp)
[2454]	168	zslope2 = zslope2 *zslope2
	169	ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) &
[2456]	170	& + zah * ( zbu * ze3wr / ( e1t(ji+ip,jj) * e2t(ji+ip,jj) ) ) * zslope2
[2450]	171	IF( ln_traldf_gdia ) THEN
[3294]	172	zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew ! fsaeit(ji+ip,jj,jk)*zslope_skew
[2454]	173	psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp * zaei_slp
[2450]	174	ENDIF
	175	END DO
	176	END DO
	177	END DO
	178	END DO
	179	END DO
	180	!
[2454]	181	DO jp = 0, 1
	182	DO kp = 0, 1
	183	DO jk = 1, jpkm1
	184	DO jj = 1, jpjm1
[2450]	185	DO ji=1,fs_jpim1
[3294]	186	ze2vr = 1._wp / e2v(ji,jj)
[2454]	187	ze3wr = 1.0_wp / fse3w(ji,jj+jp,jk+kp)
	188	zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
[3294]	189	zah = fsahtv(ji,jj,jk) ! fsaht(ji,jj+jp,jk)
[2450]	190	zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
[3294]	191	! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
	192	! (do this by adding gradient of depth)
	193	zslope2 = zslope_skew + ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * ze2vr * vmask(ji,jj,jk+kp)
[2454]	194	zslope2 = zslope2 * zslope2
	195	ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) &
	196	& + zah * ( zbv * ze3wr / ( e1t(ji,jj+jp) * e2t(ji,jj+jp) ) ) * zslope2
[2450]	197	IF( ln_traldf_gdia ) THEN
[3294]	198	zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew ! fsaeit(ji,jj+jp,jk)*zslope_skew
[2454]	199	psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp * zaei_slp
[2450]	200	ENDIF
	201	END DO
	202	END DO
	203	END DO
	204	END DO
[2205]	205	END DO
[2371]	206	!
[5147]	207	IF( iom_use("uoce_eiv") .OR. iom_use("voce_eiv") .OR. iom_use("woce_eiv") ) THEN
	208	!
	209	IF( ln_traldf_gdia .AND. cdtype == 'TRA' ) THEN
	210	CALL wrk_alloc( jpi , jpj , jpk , zw3d )
	211	DO jk=1,jpkm1
	212	zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk) ! u_eiv = -dpsix/dz
	213	END DO
	214	zw3d(:,:,jpk) = 0._wp
	215	CALL iom_put( "uoce_eiv", zw3d ) ! i-eiv current
[3294]	216
[5147]	217	DO jk=1,jpk-1
	218	zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk) ! v_eiv = -dpsiy/dz
	219	END DO
	220	zw3d(:,:,jpk) = 0._wp
	221	CALL iom_put( "voce_eiv", zw3d ) ! j-eiv current
[3294]	222
[5147]	223	DO jk=1,jpk-1
	224	DO jj = 2, jpjm1
	225	DO ji = fs_2, fs_jpim1 ! vector opt.
	226	zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2t(ji,jj) + &
	227	& (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1t(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx
	228	END DO
[3294]	229	END DO
	230	END DO
[5147]	231	zw3d(:,:,jpk) = 0._wp
	232	CALL iom_put( "woce_eiv", zw3d ) ! vert. eiv current
	233	CALL wrk_dealloc( jpi , jpj , jpk , zw3d )
	234	ENDIF
	235	!
[3294]	236	ENDIF
[2371]	237	! ! ===========
	238	DO jn = 1, kjpt ! tracer loop
	239	! ! ===========
	240	! Zero fluxes for each tracer
	241	ztfw(:,:,:) = 0._wp
	242	zftu(:,:,:) = 0._wp
	243	zftv(:,:,:) = 0._wp
[3294]	244	!
[2371]	245	DO jk = 1, jpkm1 !== before lateral T & S gradients at T-level jk ==!
	246	DO jj = 1, jpjm1
	247	DO ji = 1, fs_jpim1 ! vector opt.
	248	zdit(ji,jj,jk) = ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
	249	zdjt(ji,jj,jk) = ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
	250	END DO
[2205]	251	END DO
	252	END DO
[3294]	253	IF( ln_zps.and.l_grad_zps ) THEN ! partial steps: correction at the last level
[2371]	254	DO jj = 1, jpjm1
	255	DO ji = 1, jpim1
[3294]	256	zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
	257	zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
[2371]	258	END DO
	259	END DO
	260	ENDIF
[2205]	261
[2371]	262	!!----------------------------------------------------------------------
	263	!! II - horizontal trend (full)
	264	!!----------------------------------------------------------------------
	265	!
	266	DO jk = 1, jpkm1
	267	!
	268	! !== Vertical tracer gradient at level jk and jk+1
[3294]	269	zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
[2371]	270	!
[3294]	271	! ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1)
	272	IF( jk == 1 ) THEN ; zdkt3d(:,:,0) = zdkt3d(:,:,1)
	273	ELSE ; zdkt3d(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk)
[2371]	274	ENDIF
[2205]	275
	276
[3294]	277	IF (ln_botmix_grif) THEN
	278	DO ip = 0, 1 !== Horizontal & vertical fluxes
	279	DO kp = 0, 1
	280	DO jj = 1, jpjm1
	281	DO ji = 1, fs_jpim1
	282	ze1ur = 1._wp / e1u(ji,jj)
	283	zdxt = zdit(ji,jj,jk) * ze1ur
	284	ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
	285	zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr
	286	zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
	287	zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp)
[2205]	288
[3294]	289	zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
	290	! ln_botmix_grif is .T. don't mask zah for bottom half cells
	291	zah = fsahtu(ji,jj,jk) !*umask(ji,jj,jk+kp) !fsaht(ji+ip,jj,jk) ===>> ????
	292	zah_slp = zah * zslope_iso
	293	zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew !fsaeit(ji+ip,jj,jk)*zslope_skew
	294	zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
	295	ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr
	296	END DO
[2371]	297	END DO
	298	END DO
	299	END DO
[2205]	300
[3294]	301	DO jp = 0, 1
	302	DO kp = 0, 1
	303	DO jj = 1, jpjm1
	304	DO ji = 1, fs_jpim1
	305	ze2vr = 1._wp / e2v(ji,jj)
	306	zdyt = zdjt(ji,jj,jk) * ze2vr
	307	ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp)
	308	zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr
	309	zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
	310	zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp)
	311	zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
	312	! ln_botmix_grif is .T. don't mask zah for bottom half cells
	313	zah = fsahtv(ji,jj,jk) !*vmask(ji,jj,jk+kp) ! fsaht(ji,jj+jp,jk)
	314	zah_slp = zah * zslope_iso
	315	zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew ! fsaeit(ji,jj+jp,jk)*zslope_skew
	316	zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
	317	ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr
	318	END DO
[2371]	319	END DO
	320	END DO
	321	END DO
[3294]	322	ELSE
	323	DO ip = 0, 1 !== Horizontal & vertical fluxes
	324	DO kp = 0, 1
	325	DO jj = 1, jpjm1
	326	DO ji = 1, fs_jpim1
	327	ze1ur = 1._wp / e1u(ji,jj)
	328	zdxt = zdit(ji,jj,jk) * ze1ur
	329	ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
	330	zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr
	331	zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
	332	zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp)
[2205]	333
[3294]	334	zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
	335	! ln_botmix_grif is .F. mask zah for bottom half cells
	336	zah = fsahtu(ji,jj,jk) * umask(ji,jj,jk+kp) ! fsaht(ji+ip,jj,jk) ===>> ????
	337	zah_slp = zah * zslope_iso
	338	zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew ! fsaeit(ji+ip,jj,jk)*zslope_skew
	339	zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
	340	ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr
	341	END DO
	342	END DO
	343	END DO
	344	END DO
	345
	346	DO jp = 0, 1
	347	DO kp = 0, 1
	348	DO jj = 1, jpjm1
	349	DO ji = 1, fs_jpim1
	350	ze2vr = 1._wp / e2v(ji,jj)
	351	zdyt = zdjt(ji,jj,jk) * ze2vr
	352	ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp)
	353	zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr
	354	zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
	355	zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp)
	356	zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
	357	! ln_botmix_grif is .F. mask zah for bottom half cells
	358	zah = fsahtv(ji,jj,jk) * vmask(ji,jj,jk+kp) ! fsaht(ji,jj+jp,jk)
	359	zah_slp = zah * zslope_iso
	360	zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew ! fsaeit(ji,jj+jp,jk)*zslope_skew
	361	zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
	362	ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr
	363	END DO
	364	END DO
	365	END DO
	366	END DO
	367	END IF
	368	! !== divergence and add to the general trend ==!
[2450]	369	DO jj = 2 , jpjm1
[3294]	370	DO ji = fs_2, fs_jpim1 ! vector opt.
[2450]	371	zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	372	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * ( zftu(ji-1,jj,jk) - zftu(ji,jj,jk) &
	373	& + zftv(ji,jj-1,jk) - zftv(ji,jj,jk) )
	374	END DO
	375	END DO
	376	!
	377	END DO
	378	!
[3294]	379	DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to the general tracer trend
[2450]	380	DO jj = 2, jpjm1
[3294]	381	DO ji = fs_2, fs_jpim1 ! vector opt.
[2450]	382	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) &
	383	& / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	384	END DO
	385	END DO
	386	END DO
	387	!
[3294]	388	! ! "Poleward" diffusive heat or salt transports (T-S case only)
[7179]	389	IF( cdtype == 'TRA' .AND. ln_diaptr ) CALL dia_ptr_ohst_components( jn, 'ldf', zftv(:,:,:) )
[2205]	390
[5147]	391	IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN
	392	!
	393	IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN
	394	z2d(:,:) = 0._wp
	395	DO jk = 1, jpkm1
	396	DO jj = 2, jpjm1
	397	DO ji = fs_2, fs_jpim1 ! vector opt.
	398	z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)
	399	END DO
[2450]	400	END DO
	401	END DO
[5147]	402	z2d(:,:) = rau0_rcp * z2d(:,:)
	403	CALL lbc_lnk( z2d, 'U', -1. )
	404	CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction
	405	!
	406	z2d(:,:) = 0._wp
	407	DO jk = 1, jpkm1
	408	DO jj = 2, jpjm1
	409	DO ji = fs_2, fs_jpim1 ! vector opt.
	410	z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk)
	411	END DO
[2450]	412	END DO
	413	END DO
[5147]	414	z2d(:,:) = rau0_rcp * z2d(:,:)
	415	CALL lbc_lnk( z2d, 'V', -1. )
	416	CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction
	417	END IF
	418	!
	419	ENDIF
[2450]	420	!
	421	END DO
	422	!
[3294]	423	CALL wrk_dealloc( jpi, jpj, z2d )
	424	CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt, ztfw )
[2715]	425	!
[3294]	426	IF( nn_timing == 1 ) CALL timing_stop('tra_ldf_iso_grif')
	427	!
[2371]	428	END SUBROUTINE tra_ldf_iso_grif
	429
[2205]	430	#else
	431	!!----------------------------------------------------------------------
	432	!! default option : Dummy code NO rotation of the diffusive tensor
	433	!!----------------------------------------------------------------------
[3294]	434	REAL, PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: psix_eiv, psiy_eiv !: eiv stream function (diag only)
[2205]	435	CONTAINS
[3294]	436	SUBROUTINE tra_ldf_iso_grif( kt, kit000, cdtype, pgu, pgv, &
	437	& ptb, pta, kjpt, pahtb0 )
[2454]	438	CHARACTER(len=3) :: cdtype
[3294]	439	INTEGER :: kit000 ! first time step index
[2454]	440	REAL, DIMENSION(:,:,:) :: pgu, pgv ! tracer gradient at pstep levels
	441	REAL, DIMENSION(:,:,:,:) :: ptb, pta
	442	WRITE(,) 'tra_ldf_iso_grif: You should not have seen this print! error?', kt, cdtype, &
	443	& pgu(1,1,1), pgv(1,1,1), ptb(1,1,1,1), pta(1,1,1,1), kjpt, pahtb0
[2205]	444	END SUBROUTINE tra_ldf_iso_grif
	445	#endif
	446
	447	!!==============================================================================
	448	END MODULE traldf_iso_grif

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/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90 @ 13320

Download in other formats: