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.

zdfric.F90 in branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

source: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90 @ 7953

Last change on this file since 7953 was 7953, checked in by gm, 7 years ago
#1880 (HPC-09): add zdfphy (the ZDF manager) + remove all key_...
Property svn:keywords set to `Id`
File size: 15.2 KB

Rev	Line
[3]	1	MODULE zdfric
	2	!!======================================================================
	3	!! * MODULE zdfric *
	4	!! Ocean physics: vertical mixing coefficient compute from the local
	5	!! Richardson number dependent formulation
	6	!!======================================================================
[1601]	7	!! History : OPA ! 1987-09 (P. Andrich) Original code
	8	!! 4.0 ! 1991-11 (G. Madec)
[3294]	9	!! 7.0 ! 1996-01 (G. Madec) complete rewriting of multitasking suppression of common work arrays
	10	!! 8.0 ! 1997-06 (G. Madec) complete rewriting of zdfmix
[1601]	11	!! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module
[2528]	12	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
[3294]	13	!! 3.3.1! 2011-09 (P. Oddo) Mixed layer depth parameterization
[1601]	14	!!----------------------------------------------------------------------
[7953]	15
[3]	16	!!----------------------------------------------------------------------
[7953]	17	!! zdf_ric : update momentum and tracer Kz from the Richardson number
[3625]	18	!! zdf_ric_init : initialization, namelist read, & parameters control
[3]	19	!!----------------------------------------------------------------------
[3625]	20	USE oce ! ocean dynamics and tracers variables
	21	USE dom_oce ! ocean space and time domain variables
	22	USE zdf_oce ! ocean vertical physics
	23	USE in_out_manager ! I/O manager
	24	USE lbclnk ! ocean lateral boundary condition (or mpp link)
	25	USE lib_mpp ! MPP library
	26	USE wrk_nemo ! work arrays
	27	USE timing ! Timing
	28	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[3]	29
[6489]	30	USE eosbn2, ONLY : neos
[3294]	31
[3]	32	IMPLICIT NONE
	33	PRIVATE
	34
[2528]	35	PUBLIC zdf_ric ! called by step.F90
	36	PUBLIC zdf_ric_init ! called by opa.F90
[3]	37
[4147]	38	! !!* Namelist namzdf_ric : Richardson number dependent Kz *
	39	INTEGER :: nn_ric ! coefficient of the parameterization
	40	REAL(wp) :: rn_avmri ! maximum value of the vertical eddy viscosity
	41	REAL(wp) :: rn_alp ! coefficient of the parameterization
	42	REAL(wp) :: rn_ekmfc ! Ekman Factor Coeff
	43	REAL(wp) :: rn_mldmin ! minimum mixed layer (ML) depth
	44	REAL(wp) :: rn_mldmax ! maximum mixed layer depth
	45	REAL(wp) :: rn_wtmix ! Vertical eddy Diff. in the ML
	46	REAL(wp) :: rn_wvmix ! Vertical eddy Visc. in the ML
	47	LOGICAL :: ln_mldw ! Use or not the MLD parameters
[3]	48
[2715]	49	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tmric !: coef. for the horizontal mean at t-point
[1537]	50
[3]	51	!! * Substitutions
[6140]	52	# include "vectopt_loop_substitute.h90"
[3]	53	!!----------------------------------------------------------------------
[2715]	54	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
[2528]	55	!! $Id$
[2715]	56	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[247]	57	!!----------------------------------------------------------------------
[3]	58	CONTAINS
	59
[2715]	60	INTEGER FUNCTION zdf_ric_alloc()
	61	!!----------------------------------------------------------------------
	62	!! * FUNCTION zdf_ric_alloc *
	63	!!----------------------------------------------------------------------
	64	ALLOCATE( tmric(jpi,jpj,jpk) , STAT= zdf_ric_alloc )
	65	!
	66	IF( lk_mpp ) CALL mpp_sum ( zdf_ric_alloc )
	67	IF( zdf_ric_alloc /= 0 ) CALL ctl_warn('zdf_ric_alloc: failed to allocate arrays')
	68	END FUNCTION zdf_ric_alloc
	69
	70
[3]	71	SUBROUTINE zdf_ric( kt )
	72	!!----------------------------------------------------------------------
	73	!! * ROUTINE zdfric *
	74	!!
	75	!! ** Purpose : Compute the before eddy viscosity and diffusivity as
[3294]	76	!! a function of the local richardson number.
[3]	77	!!
	78	!! ** Method : Local richardson number dependent formulation of the
[3294]	79	!! vertical eddy viscosity and diffusivity coefficients.
[1601]	80	!! The eddy coefficients are given by:
	81	!! avm = avm0 + avmb
	82	!! avt = avm0 / (1 + rn_alp*ri)
[3294]	83	!! with ri = N^2 / dz(u)**2
	84	!! = e3w*2 rn2/[ mi( dk(ub) )+mj( dk(vb) ) ]
	85	!! avm0= rn_avmri / (1 + rn_alpri)*nn_ric
[1601]	86	!! Where ri is the before local Richardson number,
	87	!! rn_avmri is the maximum value reaches by avm and avt
	88	!! avmb and avtb are the background (or minimum) values
	89	!! and rn_alp, nn_ric are adjustable parameters.
	90	!! Typical values used are : avm0=1.e-2 m2/s, avmb=1.e-6 m2/s
[1537]	91	!! avtb=1.e-7 m2/s, rn_alp=5. and nn_ric=2.
[3]	92	!! a numerical threshold is impose on the vertical shear (1.e-20)
[3294]	93	!! As second step compute Ekman depth from wind stress forcing
	94	!! and apply namelist provided vertical coeff within this depth.
	95	!! The Ekman depth is:
	96	!! Ustar = SQRT(Taum/rho0)
	97	!! ekd= rn_ekmfc * Ustar / f0
	98	!! Large et al. (1994, eq.29) suggest rn_ekmfc=0.7; however, the derivation
	99	!! of the above equation indicates the value is somewhat arbitrary; therefore
	100	!! we allow the freedom to increase or decrease this value, if the
	101	!! Ekman depth estimate appears too shallow or too deep, respectively.
	102	!! Ekd is then limited by rn_mldmin and rn_mldmax provided in the
	103	!! namelist
[3]	104	!! N.B. the mask are required for implicit scheme, and surface
[7953]	105	!! and bottom value already set in zdfphy.F90
[3]	106	!!
[1601]	107	!! References : Pacanowski & Philander 1981, JPO, 1441-1451.
[3294]	108	!! PFJ Lermusiaux 2001.
[3]	109	!!----------------------------------------------------------------------
[3294]	110	USE phycst, ONLY: rsmall,rau0
	111	USE sbc_oce, ONLY: taum
[2715]	112	!!
[3294]	113	INTEGER, INTENT( in ) :: kt ! ocean time-step
[1537]	114	!!
[3294]	115	INTEGER :: ji, jj, jk ! dummy loop indices
	116	REAL(wp) :: zcoef, zdku, zdkv, zri, z05alp, zflageos ! temporary scalars
	117	REAL(wp) :: zrhos, zustar
	118	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, ekm_dep
[1537]	119	!!----------------------------------------------------------------------
[3294]	120	!
	121	IF( nn_timing == 1 ) CALL timing_start('zdf_ric')
	122	!
	123	CALL wrk_alloc( jpi,jpj, zwx, ekm_dep )
[3]	124	! ! ===============
	125	DO jk = 2, jpkm1 ! Horizontal slab
	126	! ! ===============
	127	! Richardson number (put in zwx(ji,jj))
	128	! -----------------
	129	DO jj = 2, jpjm1
[6140]	130	DO ji = fs_2, fs_jpim1
	131	zcoef = 0.5 / e3w_n(ji,jj,jk)
[1601]	132	! ! shear of horizontal velocity
[3]	133	zdku = zcoef * ( ub(ji-1,jj,jk-1) + ub(ji,jj,jk-1) &
[1537]	134	& -ub(ji-1,jj,jk ) - ub(ji,jj,jk ) )
[3]	135	zdkv = zcoef * ( vb(ji,jj-1,jk-1) + vb(ji,jj,jk-1) &
[1537]	136	& -vb(ji,jj-1,jk ) - vb(ji,jj,jk ) )
[1601]	137	! ! richardson number (minimum value set to zero)
[3]	138	zri = rn2(ji,jj,jk) / ( zdkuzdku + zdkvzdkv + 1.e-20 )
	139	zwx(ji,jj) = MAX( zri, 0.e0 )
	140	END DO
	141	END DO
[1601]	142	CALL lbc_lnk( zwx, 'W', 1. ) ! Boundary condition (sign unchanged)
[3]	143
	144	! Vertical eddy viscosity and diffusivity coefficients
	145	! -------------------------------------------------------
[2715]	146	z05alp = 0.5_wp * rn_alp
[1601]	147	DO jj = 1, jpjm1 ! Eddy viscosity coefficients (avm)
[6140]	148	DO ji = 1, fs_jpim1
[2715]	149	avmu(ji,jj,jk) = umask(ji,jj,jk) * rn_avmri / ( 1. + z05alp( zwx(ji+1,jj)+zwx(ji,jj) ) )*nn_ric
	150	avmv(ji,jj,jk) = vmask(ji,jj,jk) * rn_avmri / ( 1. + z05alp( zwx(ji,jj+1)+zwx(ji,jj) ) )*nn_ric
[3]	151	END DO
	152	END DO
[1601]	153	DO jj = 2, jpjm1 ! Eddy diffusivity coefficients (avt)
[6140]	154	DO ji = fs_2, fs_jpim1
[2715]	155	avt(ji,jj,jk) = tmric(ji,jj,jk) / ( 1._wp + rn_alp * zwx(ji,jj) ) &
	156	& * ( avmu(ji,jj,jk) + avmu(ji-1,jj,jk) &
	157	& + avmv(ji,jj,jk) + avmv(ji,jj-1,jk) ) &
[1537]	158	& + avtb(jk) * tmask(ji,jj,jk)
[7048]	159	END DO
	160	END DO
	161	DO jj = 2, jpjm1 ! Add the background coefficient on eddy viscosity
	162	DO ji = fs_2, fs_jpim1
[3]	163	avmu(ji,jj,jk) = avmu(ji,jj,jk) + avmb(jk) * umask(ji,jj,jk)
	164	avmv(ji,jj,jk) = avmv(ji,jj,jk) + avmb(jk) * vmask(ji,jj,jk)
	165	END DO
	166	END DO
	167	! ! ===============
	168	END DO ! End of slab
	169	! ! ===============
[1601]	170	!
[3294]	171	IF( ln_mldw ) THEN
	172
	173	! Compute Ekman depth from wind stress forcing.
	174	! -------------------------------------------------------
[6489]	175	zflageos = ( 0.5 + SIGN( 0.5, neos - 1. ) ) * rau0
[6140]	176	DO jj = 2, jpjm1
	177	DO ji = fs_2, fs_jpim1
[3294]	178	zrhos = rhop(ji,jj,1) + zflageos * ( 1. - tmask(ji,jj,1) )
	179	zustar = SQRT( taum(ji,jj) / ( zrhos + rsmall ) )
[7953]	180	ekm_dep(ji,jj) = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )
[3294]	181	ekm_dep(ji,jj) = MAX(ekm_dep(ji,jj),rn_mldmin) ! Minimun allowed
	182	ekm_dep(ji,jj) = MIN(ekm_dep(ji,jj),rn_mldmax) ! Maximum allowed
	183	END DO
	184	END DO
	185
	186	! In the first model level vertical diff/visc coeff.s
	187	! are always equal to the namelist values rn_wtmix/rn_wvmix
	188	! -------------------------------------------------------
[6140]	189	DO jj = 2, jpjm1
	190	DO ji = fs_2, fs_jpim1
[3294]	191	avmv(ji,jj,1) = MAX( avmv(ji,jj,1), rn_wvmix )
	192	avmu(ji,jj,1) = MAX( avmu(ji,jj,1), rn_wvmix )
	193	avt( ji,jj,1) = MAX( avt(ji,jj,1), rn_wtmix )
	194	END DO
	195	END DO
	196
	197	! Force the vertical mixing coef within the Ekman depth
	198	! -------------------------------------------------------
	199	DO jk = 2, jpkm1
[6140]	200	DO jj = 2, jpjm1
	201	DO ji = fs_2, fs_jpim1
[7646]	202	IF( gdept_n(ji,jj,jk) < ekm_dep(ji,jj) ) THEN
[3294]	203	avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), rn_wvmix )
	204	avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), rn_wvmix )
	205	avt( ji,jj,jk) = MAX( avt(ji,jj,jk), rn_wtmix )
	206	ENDIF
	207	END DO
	208	END DO
	209	END DO
	210
	211	DO jk = 1, jpkm1
[6140]	212	DO jj = 2, jpjm1
	213	DO ji = fs_2, fs_jpim1
[3294]	214	avmv(ji,jj,jk) = avmv(ji,jj,jk) * vmask(ji,jj,jk)
	215	avmu(ji,jj,jk) = avmu(ji,jj,jk) * umask(ji,jj,jk)
	216	avt( ji,jj,jk) = avt( ji,jj,jk) * tmask(ji,jj,jk)
	217	END DO
	218	END DO
	219	END DO
	220
	221	ENDIF
	222
[1601]	223	CALL lbc_lnk( avt , 'W', 1. ) ! Boundary conditions (unchanged sign)
	224	CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. )
	225	!
[3294]	226	CALL wrk_dealloc( jpi,jpj, zwx, ekm_dep )
[2715]	227	!
[3294]	228	IF( nn_timing == 1 ) CALL timing_stop('zdf_ric')
	229	!
[3]	230	END SUBROUTINE zdf_ric
	231
	232
	233	SUBROUTINE zdf_ric_init
	234	!!----------------------------------------------------------------------
	235	!! * ROUTINE zdfbfr_init *
	236	!!
	237	!! ** Purpose : Initialization of the vertical eddy diffusivity and
	238	!! viscosity coef. for the Richardson number dependent formulation.
	239	!!
[1601]	240	!! ** Method : Read the namzdf_ric namelist and check the parameter values
[3]	241	!!
[1601]	242	!! ** input : Namelist namzdf_ric
[3]	243	!!
	244	!! ** Action : increase by 1 the nstop flag is setting problem encounter
	245	!!----------------------------------------------------------------------
[2715]	246	INTEGER :: ji, jj, jk ! dummy loop indices
[4147]	247	INTEGER :: ios ! Local integer output status for namelist read
[1537]	248	!!
[3294]	249	NAMELIST/namzdf_ric/ rn_avmri, rn_alp , nn_ric , rn_ekmfc, &
	250	& rn_mldmin, rn_mldmax, rn_wtmix, rn_wvmix, ln_mldw
[3]	251	!!----------------------------------------------------------------------
[1601]	252	!
[4147]	253	REWIND( numnam_ref ) ! Namelist namzdf_ric in reference namelist : Vertical diffusion Kz depends on Richardson number
	254	READ ( numnam_ref, namzdf_ric, IOSTAT = ios, ERR = 901)
	255	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ric in reference namelist', lwp )
	256
	257	REWIND( numnam_cfg ) ! Namelist namzdf_ric in configuration namelist : Vertical diffusion Kz depends on Richardson number
	258	READ ( numnam_cfg, namzdf_ric, IOSTAT = ios, ERR = 902 )
	259	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ric in configuration namelist', lwp )
[4624]	260	IF(lwm) WRITE ( numond, namzdf_ric )
[1601]	261	!
	262	IF(lwp) THEN ! Control print
[3]	263	WRITE(numout,*)
[1537]	264	WRITE(numout,*) 'zdf_ric : Ri depend vertical mixing scheme'
	265	WRITE(numout,*) '~~~~~~~'
[1601]	266	WRITE(numout,*) ' Namelist namzdf_ric : set Kz(Ri) parameters'
[3294]	267	WRITE(numout,*) ' maximum vertical viscosity rn_avmri = ', rn_avmri
	268	WRITE(numout,*) ' coefficient rn_alp = ', rn_alp
	269	WRITE(numout,*) ' coefficient nn_ric = ', nn_ric
	270	WRITE(numout,*) ' Ekman Factor Coeff rn_ekmfc = ', rn_ekmfc
	271	WRITE(numout,*) ' minimum mixed layer depth rn_mldmin = ', rn_mldmin
	272	WRITE(numout,*) ' maximum mixed layer depth rn_mldmax = ', rn_mldmax
	273	WRITE(numout,*) ' Vertical eddy Diff. in the ML rn_wtmix = ', rn_wtmix
	274	WRITE(numout,*) ' Vertical eddy Visc. in the ML rn_wvmix = ', rn_wvmix
	275	WRITE(numout,*) ' Use the MLD parameterization ln_mldw = ', ln_mldw
[3]	276	ENDIF
[1601]	277	!
[2715]	278	! ! allocate zdfric arrays
	279	IF( zdf_ric_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_ric_init : unable to allocate arrays' )
	280	!
	281	DO jk = 1, jpk ! weighting mean array tmric for 4 T-points
	282	DO jj = 2, jpj ! which accounts for coastal boundary conditions
[3]	283	DO ji = 2, jpi
	284	tmric(ji,jj,jk) = tmask(ji,jj,jk) &
[1601]	285	& / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) &
	286	& + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) )
[3]	287	END DO
	288	END DO
	289	END DO
[2715]	290	tmric(:,1,:) = 0._wp
[1601]	291	!
	292	DO jk = 1, jpk ! Initialization of vertical eddy coef. to the background value
[3]	293	avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
	294	avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
	295	avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
	296	END DO
[1601]	297	!
[3]	298	END SUBROUTINE zdf_ric_init
	299
	300	!!======================================================================
	301	END MODULE zdfric

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

Download in other formats: