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.

tramle.F90 in NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/tramle.F90 @ 13516

Last change on this file since 13516 was 13516, checked in by hadcv, 4 years ago
Tiling for tra_adv
Property svn:keywords set to `Id`
File size: 17.6 KB

Rev	Line
[9531]	1	MODULE tramle
[3959]	2	!!======================================================================
[9531]	3	!! * MODULE tramle *
[3959]	4	!! Ocean tracers: Mixed Layer Eddy induced transport
	5	!!======================================================================
	6	!! History : 3.3 ! 2010-08 (G. Madec) Original code
	7	!!----------------------------------------------------------------------
	8
	9	!!----------------------------------------------------------------------
[9531]	10	!! tra_mle_trp : update the effective transport with the Mixed Layer Eddy induced transport
	11	!! tra_mle_init : initialisation of the Mixed Layer Eddy induced transport computation
[3959]	12	!!----------------------------------------------------------------------
	13	USE oce ! ocean dynamics and tracers variables
	14	USE dom_oce ! ocean space and time domain variables
	15	USE phycst ! physical constant
	16	USE zdfmxl ! mixed layer depth
[9019]	17	!
[3959]	18	USE in_out_manager ! I/O manager
	19	USE iom ! IOM library
	20	USE lib_mpp ! MPP library
[9124]	21	USE lbclnk ! lateral boundary condition / mpp link
[3959]	22
	23	IMPLICIT NONE
	24	PRIVATE
	25
[9531]	26	PUBLIC tra_mle_trp ! routine called in traadv.F90
	27	PUBLIC tra_mle_init ! routine called in traadv.F90
[3959]	28
[9531]	29	! !!* namelist namtra_mle *
	30	LOGICAL, PUBLIC :: ln_mle !: flag to activate the Mixed Layer Eddy (MLE) parameterisation
	31	INTEGER :: nn_mle ! MLE type: =0 standard Fox-Kemper ; =1 new formulation
	32	INTEGER :: nn_mld_uv ! space interpolation of MLD at u- & v-pts (0=min,1=averaged,2=max)
	33	INTEGER :: nn_conv ! =1 no MLE in case of convection ; =0 always MLE
	34	REAL(wp) :: rn_ce ! MLE coefficient
[5836]	35	! ! parameters used in nn_mle = 0 case
[9531]	36	REAL(wp) :: rn_lf ! typical scale of mixed layer front
	37	REAL(wp) :: rn_time ! time scale for mixing momentum across the mixed layer
[5836]	38	! ! parameters used in nn_mle = 1 case
[9531]	39	REAL(wp) :: rn_lat ! reference latitude for a 5 km scale of ML front
	40	REAL(wp) :: rn_rho_c_mle ! Density criterion for definition of MLD used by FK
[3959]	41
	42	REAL(wp) :: r5_21 = 5.e0 / 21.e0 ! factor used in mle streamfunction computation
[12489]	43	REAL(wp) :: rb_c ! ML buoyancy criteria = g rho_c /rho0 where rho_c is defined in zdfmld
[3959]	44	REAL(wp) :: rc_f ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_mle=1 case
	45
	46	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: rfu, rfv ! modified Coriolis parameter (f+tau) at u- & v-pts
	47	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ft ! inverse of the modified Coriolis parameter at t-pts
	48
	49	!! * Substitutions
[12377]	50	# include "do_loop_substitute.h90"
[13237]	51	# include "domzgr_substitute.h90"
[3959]	52	!!----------------------------------------------------------------------
[9598]	53	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5215]	54	!! $Id$
[10068]	55	!! Software governed by the CeCILL license (see ./LICENSE)
[3959]	56	!!----------------------------------------------------------------------
	57	CONTAINS
	58
[12377]	59	SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype, Kmm )
[3959]	60	!!----------------------------------------------------------------------
[9531]	61	!! * ROUTINE tra_mle_trp *
[3959]	62	!!
	63	!! ** Purpose : Add to the transport the Mixed Layer Eddy induced transport
	64	!!
	65	!! ** Method : The 3 components of the Mixed Layer Eddy (MLE) induced
	66	!! transport are computed as follows :
	67	!! zu_mle = dk[ zpsi_uw ]
	68	!! zv_mle = dk[ zpsi_vw ]
	69	!! zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ]
	70	!! where zpsi is the MLE streamfunction at uw and vw points (see the doc)
	71	!! and added to the input velocity :
	72	!! p.n = p.n + z._mle
	73	!!
[12377]	74	!! ** Action : - (pu,pv,pw) increased by the mle transport
[3959]	75	!! CAUTION, the transport is not updated at the last line/raw
	76	!! this may be a problem for some advection schemes
	77	!!
	78	!! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
	79	!! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
	80	!!----------------------------------------------------------------------
[13516]	81	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	82	INTEGER , INTENT(in ) :: kit000 ! first time step index
	83	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
	84	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	85	REAL(wp), DIMENSION(ST_2D(nn_hls),jpk), INTENT(inout) :: pu ! in : 3 ocean transport components
	86	REAL(wp), DIMENSION(ST_2D(nn_hls),jpk), INTENT(inout) :: pv ! out: same 3 transport components
	87	REAL(wp), DIMENSION(ST_2D(nn_hls),jpk), INTENT(inout) :: pw ! increased by the MLE induced transport
[3959]	88	!
[9019]	89	INTEGER :: ji, jj, jk ! dummy loop indices
	90	INTEGER :: ii, ij, ik, ikmax ! local integers
	91	REAL(wp) :: zcuw, zmuw, zc ! local scalar
	92	REAL(wp) :: zcvw, zmvw ! - -
[13516]	93	INTEGER , DIMENSION(ST_2D(nn_hls)) :: inml_mle
	94	REAL(wp), DIMENSION(ST_2D(nn_hls)) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_MH
	95	REAL(wp), DIMENSION(ST_2D(nn_hls),jpk) :: zpsi_uw, zpsi_vw
	96	! TEMP: These changes not necessary if using XIOS (subdomain support)
	97	REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: zLf_NH
	98	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zpsiu_mle, zpsiv_mle
[3959]	99	!!----------------------------------------------------------------------
[5836]	100	!
[3995]	101	! !== MLD used for MLE ==!
	102	! ! compute from the 10m density to deal with the diurnal cycle
[13516]	103	DO_2D( 1, 1, 1, 1 )
	104	inml_mle(ji,jj) = mbkt(ji,jj) + 1 ! init. to number of ocean w-level (T-level + 1)
	105	END_2D
[10105]	106	IF ( nla10 > 0 ) THEN ! avoid case where first level is thicker than 10m
[13295]	107	DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )
[12377]	108	IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle ) inml_mle(ji,jj) = jk ! Mixed layer
	109	END_3D
[10105]	110	ENDIF
[4325]	111	ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 ) ! max level of the computation
[3959]	112	!
	113	!
[3995]	114	zmld(:,:) = 0._wp !== Horizontal shape of the MLE ==!
	115	zbm (:,:) = 0._wp
	116	zn2 (:,:) = 0._wp
[13295]	117	DO_3D( 1, 1, 1, 1, 1, ikmax )
[12377]	118	zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points
	119	zmld(ji,jj) = zmld(ji,jj) + zc
[12489]	120	zbm (ji,jj) = zbm (ji,jj) + zc * (rho0 - rhop(ji,jj,jk) ) * r1_rho0
[12377]	121	zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp
	122	END_3D
[3959]	123
[3995]	124	SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts
[3959]	125	CASE ( 0 ) != min of the 2 neighbour MLDs
[13295]	126	DO_2D( 1, 0, 1, 0 )
[12377]	127	zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) )
	128	zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) )
	129	END_2D
[3959]	130	CASE ( 1 ) != average of the 2 neighbour MLDs
[13295]	131	DO_2D( 1, 0, 1, 0 )
[12377]	132	zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp
	133	zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp
	134	END_2D
[3959]	135	CASE ( 2 ) != max of the 2 neighbour MLDs
[13295]	136	DO_2D( 1, 0, 1, 0 )
[12377]	137	zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) )
	138	zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) )
	139	END_2D
[3959]	140	END SELECT
[3995]	141	! ! convert density into buoyancy
[13516]	142	DO_2D( 1, 1, 1, 1 )
	143	zbm(ji,jj) = + grav * zbm(ji,jj) / MAX( e3t(ji,jj,1,Kmm), zmld(ji,jj) )
	144	END_2D
[3959]	145	!
	146	!
[3995]	147	! !== Magnitude of the MLE stream function ==!
	148	!
[3959]	149	! di[bm] Ds
	150	! Psi = Ce H^2 ---------------- e2u mu(z) where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 )
	151	! e1u Lf fu and the e2u for the "transport"
	152	! (not *e3u as divided by e3u at the end)
	153	!
	154	IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation
[13295]	155	DO_2D( 1, 0, 1, 0 )
[12377]	156	zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) &
	157	& * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) &
	158	& / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) )
	159	!
	160	zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) &
	161	& * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) &
	162	& / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) )
	163	END_2D
[3959]	164	!
	165	ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
[13295]	166	DO_2D( 1, 0, 1, 0 )
[12377]	167	zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) &
	168	& * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )
	169	!
	170	zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) &
	171	& * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )
	172	END_2D
[3959]	173	ENDIF
	174	!
	175	IF( nn_conv == 1 ) THEN ! No MLE in case of convection
[13295]	176	DO_2D( 1, 0, 1, 0 )
[12377]	177	IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp ) zpsim_u(ji,jj) = 0._wp
	178	IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp ) zpsim_v(ji,jj) = 0._wp
	179	END_2D
[3959]	180	ENDIF
	181	!
[3995]	182	! !== structure function value at uw- and vw-points ==!
[13295]	183	DO_2D( 1, 0, 1, 0 )
[12377]	184	zhu(ji,jj) = 1._wp / zhu(ji,jj) ! hu --> 1/hu
	185	zhv(ji,jj) = 1._wp / zhv(ji,jj)
	186	END_2D
[4835]	187	!
[3959]	188	zpsi_uw(:,:,:) = 0._wp
	189	zpsi_vw(:,:,:) = 0._wp
	190	!
[13295]	191	DO_3D( 1, 0, 1, 0, 2, ikmax )
[12377]	192	zcuw = 1._wp - ( gdepw(ji+1,jj,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhu(ji,jj)
	193	zcvw = 1._wp - ( gdepw(ji,jj+1,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhv(ji,jj)
	194	zcuw = zcuw * zcuw
	195	zcvw = zcvw * zcvw
	196	zmuw = MAX( 0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw ) )
	197	zmvw = MAX( 0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw ) )
	198	!
	199	zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk)
	200	zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk)
	201	END_3D
[3995]	202	!
	203	! !== transport increased by the MLE induced transport ==!
[3959]	204	DO jk = 1, ikmax
[13295]	205	DO_2D( 1, 0, 1, 0 )
[12377]	206	pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
	207	pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
	208	END_2D
[13295]	209	DO_2D( 0, 0, 0, 0 )
[12377]	210	pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk) &
	211	& + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) )
	212	END_2D
[3959]	213	END DO
	214
[13516]	215	! TEMP: These changes not necessary if using XIOS (subdomain support)
[3995]	216	IF( cdtype == 'TRA') THEN !== outputs ==!
[13516]	217	IF( kt == nit000 .AND. (ntile == 0 .OR. ntile == 1) ) THEN ! Do only on the first tile and timestep
	218	ALLOCATE( zLf_NH(jpi,jpj), zpsiu_mle(jpi,jpj,jpk), zpsiv_mle(jpi,jpj,jpk) )
	219	ENDIF
[3959]	220	!
[13516]	221	DO_2D( 1, 1, 1, 1 )
	222	zLf_NH(ji,jj) = SQRT( rb_c * zmld(ji,jj) ) * r1_ft(ji,jj) ! Lf = N H / f
	223	END_2D
[3959]	224	!
[4188]	225	! divide by cross distance to give streamfunction with dimensions m^2/s
[13516]	226	DO_3D( 1, 1, 1, 1, 1, ikmax+1 )
	227	zpsiu_mle(ji,jj,jk) = zpsi_uw(ji,jj,jk) * r1_e2u(ji,jj)
	228	zpsiv_mle(ji,jj,jk) = zpsi_vw(ji,jj,jk) * r1_e1v(ji,jj)
	229	END_3D
	230
	231	IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile
	232	CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f
	233	CALL iom_put( "psiu_mle", zpsiu_mle ) ! i-mle streamfunction
	234	CALL iom_put( "psiv_mle", zpsiv_mle ) ! j-mle streamfunction
	235	ENDIF
[3959]	236	ENDIF
	237	!
[9531]	238	END SUBROUTINE tra_mle_trp
[3959]	239
	240
[9531]	241	SUBROUTINE tra_mle_init
[3959]	242	!!---------------------------------------------------------------------
[9531]	243	!! * ROUTINE tra_mle_init *
[3959]	244	!!
	245	!! ** Purpose : Control the consistency between namelist options for
	246	!! tracer advection schemes and set nadv
	247	!!----------------------------------------------------------------------
	248	INTEGER :: ji, jj, jk ! dummy loop indices
	249	INTEGER :: ierr
[4245]	250	INTEGER :: ios ! Local integer output status for namelist read
[3959]	251	REAL(wp) :: z1_t2, zfu, zfv ! - -
	252	!
[9531]	253	NAMELIST/namtra_mle/ ln_mle , nn_mle, rn_ce, rn_lf, rn_time, rn_lat, nn_mld_uv, nn_conv, rn_rho_c_mle
[3959]	254	!!----------------------------------------------------------------------
	255
[9531]	256	READ ( numnam_ref, namtra_mle, IOSTAT = ios, ERR = 901)
[11536]	257	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_mle in reference namelist' )
[4245]	258
[9531]	259	READ ( numnam_cfg, namtra_mle, IOSTAT = ios, ERR = 902 )
[11536]	260	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_mle in configuration namelist' )
[9531]	261	IF(lwm) WRITE ( numond, namtra_mle )
[4245]	262
[3959]	263	IF(lwp) THEN ! Namelist print
	264	WRITE(numout,*)
[9531]	265	WRITE(numout,*) 'tra_mle_init : mixed layer eddy (MLE) advection acting on tracers'
	266	WRITE(numout,*) '~~~~~~~~~~~~~'
	267	WRITE(numout,*) ' Namelist namtra_mle : mixed layer eddy advection applied on tracers'
[9168]	268	WRITE(numout,*) ' use mixed layer eddy (MLE, i.e. Fox-Kemper param) (T/F) ln_mle = ', ln_mle
	269	WRITE(numout,*) ' MLE type: =0 standard Fox-Kemper ; =1 new formulation nn_mle = ', nn_mle
	270	WRITE(numout,*) ' magnitude of the MLE (typical value: 0.06 to 0.08) rn_ce = ', rn_ce
	271	WRITE(numout,*) ' scale of ML front (ML radius of deformation) (rn_mle=0) rn_lf = ', rn_lf, 'm'
	272	WRITE(numout,*) ' maximum time scale of MLE (rn_mle=0) rn_time = ', rn_time, 's'
	273	WRITE(numout,*) ' reference latitude (degrees) of MLE coef. (rn_mle=1) rn_lat = ', rn_lat, 'deg'
	274	WRITE(numout,*) ' space interp. of MLD at u-(v-)pts (0=min,1=averaged,2=max) nn_mld_uv = ', nn_mld_uv
	275	WRITE(numout,*) ' =1 no MLE in case of convection ; =0 always MLE nn_conv = ', nn_conv
	276	WRITE(numout,*) ' Density difference used to define ML for FK rn_rho_c_mle = ', rn_rho_c_mle
[3959]	277	ENDIF
	278	!
	279	IF(lwp) THEN
	280	WRITE(numout,*)
	281	IF( ln_mle ) THEN
[9190]	282	WRITE(numout,*) ' ==>>> Mixed Layer Eddy induced transport added to tracer advection'
[7646]	283	IF( nn_mle == 0 ) WRITE(numout,*) ' Fox-Kemper et al 2010 formulation'
	284	IF( nn_mle == 1 ) WRITE(numout,*) ' New formulation'
[3959]	285	ELSE
[9190]	286	WRITE(numout,*) ' ==>>> Mixed Layer Eddy parametrisation NOT used'
[3959]	287	ENDIF
	288	ENDIF
	289	!
	290	IF( ln_mle ) THEN ! MLE initialisation
	291	!
[12489]	292	rb_c = grav * rn_rho_c_mle /rho0 ! Mixed Layer buoyancy criteria
[3959]	293	IF(lwp) WRITE(numout,*)
	294	IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 '
	295	IF(lwp) WRITE(numout,*) ' associated ML density criteria defined in zdfmxl = ', rho_c, 'kg/m3'
	296	!
	297	IF( nn_mle == 0 ) THEN ! MLE array allocation & initialisation
	298	ALLOCATE( rfu(jpi,jpj) , rfv(jpi,jpj) , STAT= ierr )
	299	IF( ierr /= 0 ) CALL ctl_stop( 'tra_adv_mle_init: failed to allocate arrays' )
	300	z1_t2 = 1._wp / ( rn_time * rn_time )
[13295]	301	DO_2D( 0, 1, 0, 1 )
[12377]	302	zfu = ( ff_f(ji,jj) + ff_f(ji,jj-1) ) * 0.5_wp
	303	zfv = ( ff_f(ji,jj) + ff_f(ji-1,jj) ) * 0.5_wp
	304	rfu(ji,jj) = SQRT( zfu * zfu + z1_t2 )
	305	rfv(ji,jj) = SQRT( zfv * zfv + z1_t2 )
	306	END_2D
[13226]	307	CALL lbc_lnk_multi( 'tramle', rfu, 'U', 1.0_wp , rfv, 'V', 1.0_wp )
[3959]	308	!
	309	ELSEIF( nn_mle == 1 ) THEN ! MLE array allocation & initialisation
	310	rc_f = rn_ce / ( 5.e3_wp * 2._wp * omega * SIN( rad * rn_lat ) )
	311	!
	312	ENDIF
	313	!
	314	! ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_mle case)
	315	ALLOCATE( r1_ft(jpi,jpj) , STAT= ierr )
	316	IF( ierr /= 0 ) CALL ctl_stop( 'tra_adv_mle_init: failed to allocate r1_ft array' )
	317	!
	318	z1_t2 = 1._wp / ( rn_time * rn_time )
[7753]	319	r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 )
[3959]	320	!
	321	ENDIF
	322	!
[9531]	323	END SUBROUTINE tra_mle_init
[3959]	324
	325	!!==============================================================================
[9531]	326	END MODULE tramle

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

Download in other formats: