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.

dynldf_iso.F90 in NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DYN/dynldf_iso.F90 @ 14644

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

Rev	Line
[3]	1	MODULE dynldf_iso
	2	!!======================================================================
	3	!! * MODULE dynldf_iso *
[5836]	4	!! Ocean dynamics: lateral viscosity trend (rotated laplacian operator)
[3]	5	!!======================================================================
[2715]	6	!! History : OPA ! 97-07 (G. Madec) Original code
	7	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
	8	!! - ! 2004-08 (C. Talandier) New trends organization
	9	!! 2.0 ! 2005-11 (G. Madec) s-coordinate: horizontal diffusion
[5836]	10	!! 3.7 ! 2014-01 (F. Lemarie, G. Madec) restructuration/simplification of ahm specification,
	11	!! ! add velocity dependent coefficient and optional read in file
[2715]	12	!!----------------------------------------------------------------------
[5836]	13
[3]	14	!!----------------------------------------------------------------------
	15	!! dyn_ldf_iso : update the momentum trend with the horizontal part
	16	!! of the lateral diffusion using isopycnal or horizon-
	17	!! tal s-coordinate laplacian operator.
	18	!!----------------------------------------------------------------------
	19	USE oce ! ocean dynamics and tracers
	20	USE dom_oce ! ocean space and time domain
[5836]	21	USE ldfdyn ! lateral diffusion: eddy viscosity coef.
	22	USE ldftra ! lateral physics: eddy diffusivity
[3]	23	USE zdf_oce ! ocean vertical physics
	24	USE ldfslp ! iso-neutral slopes
[4990]	25	!
[3]	26	USE in_out_manager ! I/O manager
[2715]	27	USE lib_mpp ! MPP library
[5836]	28	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[258]	29	USE prtctl ! Print control
[3]	30
	31	IMPLICIT NONE
	32	PRIVATE
	33
[2715]	34	PUBLIC dyn_ldf_iso ! called by step.F90
	35	PUBLIC dyn_ldf_iso_alloc ! called by nemogcm.F90
[3]	36
[9019]	37	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akzu, akzv !: vertical component of rotated lateral viscosity
	38
[2715]	39	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u ! 2D workspace (dyn_ldf_iso)
	40	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v ! - -
	41
[3]	42	!! * Substitutions
[12377]	43	# include "do_loop_substitute.h90"
[13237]	44	# include "domzgr_substitute.h90"
[14219]	45	# include "single_precision_substitute.h90"
[3]	46	!!----------------------------------------------------------------------
[9598]	47	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1156]	48	!! $Id$
[10068]	49	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	50	!!----------------------------------------------------------------------
	51	CONTAINS
	52
[2715]	53	INTEGER FUNCTION dyn_ldf_iso_alloc()
	54	!!----------------------------------------------------------------------
	55	!! * ROUTINE dyn_ldf_iso_alloc *
	56	!!----------------------------------------------------------------------
[9019]	57	ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) , &
	58	& akzv(jpi,jpj,jpk) , zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc )
[2715]	59	!
	60	IF( dyn_ldf_iso_alloc /= 0 ) CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.')
	61	END FUNCTION dyn_ldf_iso_alloc
	62
	63
[12377]	64	SUBROUTINE dyn_ldf_iso( kt, Kbb, Kmm, puu, pvv, Krhs )
[3]	65	!!----------------------------------------------------------------------
	66	!! * ROUTINE dyn_ldf_iso *
	67	!!
[455]	68	!! ** Purpose : Compute the before trend of the rotated laplacian
	69	!! operator of lateral momentum diffusion except the diagonal
	70	!! vertical term that will be computed in dynzdf module. Add it
	71	!! to the general trend of momentum equation.
[3]	72	!!
	73	!! ** Method :
[455]	74	!! The momentum lateral diffusive trend is provided by a 2nd
	75	!! order operator rotated along neutral or geopotential surfaces
	76	!! (in s-coordinates).
[3]	77	!! It is computed using before fields (forward in time) and isopyc-
[455]	78	!! nal or geopotential slopes computed in routine ldfslp.
[3]	79	!! Here, u and v components are considered as 2 independent scalar
	80	!! fields. Therefore, the property of splitting divergent and rota-
	81	!! tional part of the flow of the standard, z-coordinate laplacian
	82	!! momentum diffusion is lost.
	83	!! horizontal fluxes associated with the rotated lateral mixing:
	84	!! u-component:
[12377]	85	!! ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t di[ uu ]
	86	!! - ahmt e2t * mi-1(uslp) dk[ mi(mk(uu)) ]
	87	!! zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f dj[ uu ]
	88	!! - ahmf e1f * mi(vslp) dk[ mj(mk(uu)) ]
[3]	89	!! v-component:
[12377]	90	!! zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t di[ vv ]
	91	!! - ahmf e2t * mj(uslp) dk[ mi(mk(vv)) ]
	92	!! zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f dj[ vv ]
	93	!! - ahmt e1f * mj-1(vslp) dk[ mj(mk(vv)) ]
[3]	94	!! take the horizontal divergence of the fluxes:
	95	!! diffu = 1/(e1ue2ue3u) { di [ ziut ] + dj-1[ zjuf ] }
	96	!! diffv = 1/(e1ve2ve3v) { di-1[ zivf ] + dj [ zjvt ] }
[12377]	97	!! Add this trend to the general trend (uu(rhs),vv(rhs)):
	98	!! uu(rhs) = uu(rhs) + diffu
[455]	99	!! CAUTION: here the isopycnal part is with a coeff. of aht. This
	100	!! should be modified for applications others than orca_r2 (!!bug)
[3]	101	!!
	102	!! ** Action :
[12377]	103	!! -(puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the before geopotential harmonic mixing trend
[9019]	104	!! -(akzu,akzv) to accompt for the diagonal vertical component
	105	!! of the rotated operator in dynzdf module
[3]	106	!!----------------------------------------------------------------------
[12377]	107	INTEGER , INTENT( in ) :: kt ! ocean time-step index
	108	INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs ! ocean time level indices
[14219]	109	REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
[2715]	110	!
	111	INTEGER :: ji, jj, jk ! dummy loop indices
[9019]	112	REAL(wp) :: zabe1, zmskt, zmkt, zuav, zuwslpi, zuwslpj ! local scalars
	113	REAL(wp) :: zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj ! - -
[9490]	114	REAL(wp) :: zcof0, zcof1, zcof2, zcof3, zcof4, zaht_0 ! - -
[9019]	115	REAL(wp), DIMENSION(jpi,jpj) :: ziut, zivf, zdku, zdk1u ! 2D workspace
	116	REAL(wp), DIMENSION(jpi,jpj) :: zjuf, zjvt, zdkv, zdk1v ! - -
[2715]	117	!!----------------------------------------------------------------------
[3294]	118	!
[3]	119	IF( kt == nit000 ) THEN
	120	IF(lwp) WRITE(numout,*)
	121	IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
	122	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate horizontal diffusive operator'
[2715]	123	! ! allocate dyn_ldf_bilap arrays
	124	IF( dyn_ldf_iso_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
[3]	125	ENDIF
[216]	126
[5836]	127	!!gm bug is dyn_ldf_iso called before tra_ldf_iso .... <<<<<===== TO BE CHECKED
	128	! s-coordinate: Iso-level diffusion on momentum but not on tracer
[455]	129	IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
[2715]	130	!
[13497]	131	DO_3D( 0, 0, 0, 0, 1, jpk ) ! set the slopes of iso-level
[12377]	132	uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
	133	vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
	134	wslpi(ji,jj,jk) = - ( gdepw(ji+1,jj,jk,Kbb) - gdepw(ji-1,jj,jk,Kbb) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5
	135	wslpj(ji,jj,jk) = - ( gdepw(ji,jj+1,jk,Kbb) - gdepw(ji,jj-1,jk,Kbb) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5
	136	END_3D
[455]	137	! Lateral boundary conditions on the slopes
[14644]	138	CALL lbc_lnk( 'dynldf_iso', uslp , 'U', -1.0_wp, vslp , 'V', -1.0_wp, wslpi, 'W', -1.0_wp, wslpj, 'W', -1.0_wp )
[9019]	139	!
	140	ENDIF
[9490]	141
	142	zaht_0 = 0.5_wp * rn_Ud * rn_Ld ! aht_0 from namtra_ldf = zaht_max
	143
[3]	144	! ! ===============
	145	DO jk = 1, jpkm1 ! Horizontal slab
	146	! ! ===============
	147
	148	! Vertical u- and v-shears at level jk and jk+1
	149	! ---------------------------------------------
	150	! surface boundary condition: zdku(jk=1)=zdku(jk=2)
	151	! zdkv(jk=1)=zdkv(jk=2)
	152
[12377]	153	zdk1u(:,:) = ( puu(:,:,jk,Kbb) -puu(:,:,jk+1,Kbb) ) * umask(:,:,jk+1)
	154	zdk1v(:,:) = ( pvv(:,:,jk,Kbb) -pvv(:,:,jk+1,Kbb) ) * vmask(:,:,jk+1)
[3]	155
	156	IF( jk == 1 ) THEN
	157	zdku(:,:) = zdk1u(:,:)
	158	zdkv(:,:) = zdk1v(:,:)
	159	ELSE
[12377]	160	zdku(:,:) = ( puu(:,:,jk-1,Kbb) - puu(:,:,jk,Kbb) ) * umask(:,:,jk)
	161	zdkv(:,:) = ( pvv(:,:,jk-1,Kbb) - pvv(:,:,jk,Kbb) ) * vmask(:,:,jk)
[3]	162	ENDIF
	163
	164	! -----f-----
	165	! Horizontal fluxes on U \|
	166	! --------------------=== t u t
	167	! \|
	168	! i-flux at t-point -----f-----
	169
[455]	170	IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
[14644]	171	DO_2D( 0, 1, 0, 0 )
[13237]	172	zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) &
	173	& * MIN( e3u(ji ,jj,jk,Kmm), &
	174	& e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj)
[3]	175
[12377]	176	zmskt = 1._wp / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &
	177	& + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ) , 1._wp )
[3]	178
[12377]	179	zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
	180
	181	ziut(ji,jj) = ( zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) ) &
	182	& + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
	183	& +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)
	184	END_2D
[455]	185	ELSE ! other coordinate system (zco or sco) : e3t
[14644]	186	DO_2D( 0, 1, 0, 0 )
[13237]	187	zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) &
	188	& * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj)
[3]	189
[12377]	190	zmskt = 1._wp / MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk+1) &
	191	& + umask(ji-1,jj,jk+1) + umask(ji,jj,jk ) , 1._wp )
[455]	192
[12377]	193	zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
[455]	194
[12377]	195	ziut(ji,jj) = ( zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) ) &
	196	& + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
	197	& +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)
	198	END_2D
[455]	199	ENDIF
[3]	200
	201	! j-flux at f-point
[13295]	202	DO_2D( 1, 0, 1, 0 )
[13237]	203	zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) &
	204	& * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj)
[3]	205
[12377]	206	zmskf = 1._wp / MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) &
	207	& + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ) , 1._wp )
[3]	208
[12377]	209	zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
[3]	210
[12377]	211	zjuf(ji,jj) = ( zabe2 * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) ) &
	212	& + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) &
	213	& +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) * fmask(ji,jj,jk)
	214	END_2D
[3]	215
	216	! \| t \|
	217	! Horizontal fluxes on V \| \|
	218	! --------------------=== f---v---f
	219	! \| \|
	220	! i-flux at f-point \| t \|
	221
[14644]	222	DO_2D( 1, 0, 0, 0 )
[13237]	223	zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) &
	224	& * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj)
[3]	225
[12377]	226	zmskf = 1._wp / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &
	227	& + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ) , 1._wp )
[3]	228
[12377]	229	zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
[3]	230
[12377]	231	zivf(ji,jj) = ( zabe1 * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) ) &
	232	& + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj) &
	233	& + zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk)
	234	END_2D
[3]	235
	236	! j-flux at t-point
[455]	237	IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
[14644]	238	DO_2D( 1, 0, 0, 1 )
[13237]	239	zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) &
	240	& * MIN( e3v(ji,jj ,jk,Kmm), &
	241	& e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj)
[3]	242
[12377]	243	zmskt = 1._wp / MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	244	& + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ) , 1._wp )
[3]	245
[12377]	246	zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
[3]	247
[12377]	248	zjvt(ji,jj) = ( zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) ) &
	249	& + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &
	250	& +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)
	251	END_2D
[455]	252	ELSE ! other coordinate system (zco or sco) : e3t
[14644]	253	DO_2D( 1, 0, 0, 1 )
[13237]	254	zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) &
	255	& * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj)
[3]	256
[12377]	257	zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	258	& + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. )
[3]	259
[12377]	260	zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
[455]	261
[12377]	262	zjvt(ji,jj) = ( zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) ) &
	263	& + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &
	264	& +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)
	265	END_2D
[455]	266	ENDIF
	267
	268
[3]	269	! Second derivative (divergence) and add to the general trend
	270	! -----------------------------------------------------------
[13497]	271	DO_2D( 0, 0, 0, 0 ) !!gm Question vectop possible??? !!bug
[12377]	272	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( ziut(ji+1,jj) - ziut(ji,jj ) &
[13237]	273	& + zjuf(ji ,jj) - zjuf(ji,jj-1) ) * r1_e1e2u(ji,jj) &
	274	& / e3u(ji,jj,jk,Kmm)
[12377]	275	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zivf(ji,jj ) - zivf(ji-1,jj) &
[13237]	276	& + zjvt(ji,jj+1) - zjvt(ji,jj ) ) * r1_e1e2v(ji,jj) &
	277	& / e3v(ji,jj,jk,Kmm)
[12377]	278	END_2D
[3]	279	! ! ===============
	280	END DO ! End of slab
	281	! ! ===============
[216]	282
[455]	283	! print sum trends (used for debugging)
[14219]	284	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=CASTWP(puu(:,:,:,Krhs)), clinfo1=' ldfh - Ua: ', mask1=umask, &
	285	& tab3d_2=CASTWP(pvv(:,:,:,Krhs)), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[216]	286
	287
[455]	288	! ! ===============
	289	DO jj = 2, jpjm1 ! Vertical slab
	290	! ! ===============
	291
	292
	293	! I. vertical trends associated with the lateral mixing
	294	! =====================================================
	295	! (excluding the vertical flux proportional to dk[t]
	296
	297
	298	! I.1 horizontal momentum gradient
	299	! --------------------------------
	300
	301	DO jk = 1, jpk
	302	DO ji = 2, jpi
	303	! i-gradient of u at jj
[12377]	304	zdiu (ji,jk) = tmask(ji,jj ,jk) * ( puu(ji,jj ,jk,Kbb) - puu(ji-1,jj ,jk,Kbb) )
[455]	305	! j-gradient of u and v at jj
[12377]	306	zdju (ji,jk) = fmask(ji,jj ,jk) * ( puu(ji,jj+1,jk,Kbb) - puu(ji ,jj ,jk,Kbb) )
	307	zdjv (ji,jk) = tmask(ji,jj ,jk) * ( pvv(ji,jj ,jk,Kbb) - pvv(ji ,jj-1,jk,Kbb) )
[455]	308	! j-gradient of u and v at jj+1
[12377]	309	zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( puu(ji,jj ,jk,Kbb) - puu(ji ,jj-1,jk,Kbb) )
	310	zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pvv(ji,jj+1,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) )
[455]	311	END DO
	312	END DO
	313	DO jk = 1, jpk
	314	DO ji = 1, jpim1
	315	! i-gradient of v at jj
[12377]	316	zdiv (ji,jk) = fmask(ji,jj ,jk) * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) )
[455]	317	END DO
	318	END DO
	319
	320
	321	! I.2 Vertical fluxes
	322	! -------------------
	323
	324	! Surface and bottom vertical fluxes set to zero
	325	DO ji = 1, jpi
	326	zfuw(ji, 1 ) = 0.e0
	327	zfvw(ji, 1 ) = 0.e0
	328	zfuw(ji,jpk) = 0.e0
	329	zfvw(ji,jpk) = 0.e0
	330	END DO
	331
	332	! interior (2=<jk=<jpk-1) on U field
	333	DO jk = 2, jpkm1
	334	DO ji = 2, jpim1
[9490]	335	zcof0 = 0.5_wp * zaht_0 * umask(ji,jj,jk)
[5836]	336	!
[9019]	337	zuwslpi = zcof0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
	338	zuwslpj = zcof0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
[5836]	339	!
[9019]	340	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) &
	341	+ tmask(ji,jj,jk )+tmask(ji+1,jj,jk ) , 1. )
	342	zmkf = 1./MAX( fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1) &
	343	+ fmask(ji,jj-1,jk ) + fmask(ji,jj,jk ) , 1. )
[455]	344
[9019]	345	zcof3 = - e2u(ji,jj) * zmkt * zuwslpi
	346	zcof4 = - e1u(ji,jj) * zmkf * zuwslpj
[455]	347	! vertical flux on u field
[9019]	348	zfuw(ji,jk) = zcof3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) &
	349	& + zdiu (ji,jk ) + zdiu (ji+1,jk ) ) &
	350	& + zcof4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) &
	351	& + zdj1u(ji,jk ) + zdju (ji ,jk ) )
	352	! vertical mixing coefficient (akzu)
[9490]	353	! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
	354	akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / zaht_0
[455]	355	END DO
	356	END DO
	357
	358	! interior (2=<jk=<jpk-1) on V field
	359	DO jk = 2, jpkm1
	360	DO ji = 2, jpim1
[9490]	361	zcof0 = 0.5_wp * zaht_0 * vmask(ji,jj,jk)
[9019]	362	!
	363	zvwslpi = zcof0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
	364	zvwslpj = zcof0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
	365	!
	366	zmkf = 1./MAX( fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1) &
	367	& + fmask(ji-1,jj,jk )+fmask(ji,jj,jk ) , 1. )
	368	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1) &
	369	& + tmask(ji,jj,jk )+tmask(ji,jj+1,jk ) , 1. )
[455]	370
[9019]	371	zcof3 = - e2v(ji,jj) * zmkf * zvwslpi
	372	zcof4 = - e1v(ji,jj) * zmkt * zvwslpj
[455]	373	! vertical flux on v field
[9019]	374	zfvw(ji,jk) = zcof3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) &
	375	& + zdiv (ji,jk ) + zdiv (ji-1,jk ) ) &
	376	& + zcof4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) &
	377	& + zdjv (ji,jk ) + zdj1v(ji ,jk ) )
	378	! vertical mixing coefficient (akzv)
[9490]	379	! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
	380	akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / zaht_0
[455]	381	END DO
	382	END DO
	383
	384
	385	! I.3 Divergence of vertical fluxes added to the general tracer trend
	386	! -------------------------------------------------------------------
	387	DO jk = 1, jpkm1
	388	DO ji = 2, jpim1
[13237]	389	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) &
	390	& / e3u(ji,jj,jk,Kmm)
	391	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) &
	392	& / e3v(ji,jj,jk,Kmm)
[455]	393	END DO
	394	END DO
	395	! ! ===============
	396	END DO ! End of slab
	397	! ! ===============
[3]	398	END SUBROUTINE dyn_ldf_iso
	399
	400	!!======================================================================
	401	END MODULE dynldf_iso
[14219]	402

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

Download in other formats: