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 trunk/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90 @ 5021

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

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

Download in other formats: