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.

p4zsink.F90 in branches/CMIP5_IPSL/NEMO/TOP_SRC/PISCES – NEMO

Context Navigation

source: branches/CMIP5_IPSL/NEMO/TOP_SRC/PISCES/p4zsink.F90 @ 1808

Last change on this file since 1808 was 1808, checked in by cetlod, 14 years ago
update TOP_SRC component on CMIP5_IPSL branch to take into account bugfixes (ie vertical diffusion routines), re-organization of restart part, damping of passive tracers on closed seas for PISCES
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 30.8 KB

Rev	Line
[935]	1	MODULE p4zsink
	2	!!======================================================================
	3	!! * MODULE p4zsink *
	4	!! TOP : PISCES Compute vertical flux of particulate matter due to gravitational sinking
	5	!!======================================================================
	6	!! History : 1.0 ! 2004 (O. Aumont) Original code
	7	!! 2.0 ! 2007-12 (C. Ethe, G. Madec) F90
	8	#if defined key_pisces
	9	!!----------------------------------------------------------------------
	10	!! p4z_sink : Compute vertical flux of particulate matter due to gravitational sinking
	11	!!----------------------------------------------------------------------
	12	USE trc
	13	USE oce_trc !
[1073]	14	USE sms_pisces
[935]	15	USE prtctl_trc
[1457]	16	USE iom
[935]	17
	18	IMPLICIT NONE
	19	PRIVATE
	20
	21	PUBLIC p4z_sink ! called in p4zbio.F90
	22
	23	!! * Shared module variables
	24	REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: & !:
	25	wsbio3, wsbio4, & !: POC and GOC sinking speeds
	26	wscal !: Calcite and BSi sinking speeds
	27
	28	!! * Module variables
[1180]	29	REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: & !:
[935]	30	sinking, sinking2, & !: POC sinking fluxes (different meanings depending on the parameterization
	31	sinkcal, sinksil, & !: CaCO3 and BSi sinking fluxes
	32	sinkfer !: Small BFe sinking flux
	33
[1264]	34	REAL(wp) :: &
	35	xstep , xstep2 !: Time step duration for biology
	36
[1457]	37	INTEGER :: &
	38	iksed = 10 !
	39
[935]	40	#if defined key_kriest
	41	REAL(wp) :: &
	42	xkr_sfact = 250. , & !: Sinking factor
	43	xkr_stick = 0.2 , & !: Stickiness
	44	xkr_nnano = 2.337 , & !: Nbr of cell in nano size class
	45	xkr_ndiat = 3.718 , & !: Nbr of cell in diatoms size class
	46	xkr_nmeso = 7.147 , & !: Nbr of cell in mesozoo size class
	47	xkr_naggr = 9.877 !: Nbr of cell in aggregates size class
	48
	49	REAL(wp) :: &
	50	xkr_frac
	51
	52	REAL(wp), PUBLIC :: &
	53	xkr_dnano , & !: Size of particles in nano pool
	54	xkr_ddiat , & !: Size of particles in diatoms pool
	55	xkr_dmeso , & !: Size of particles in mesozoo pool
	56	xkr_daggr , & !: Size of particles in aggregates pool
	57	xkr_wsbio_min , & !: min vertical particle speed
	58	xkr_wsbio_max !: max vertical particle speed
	59
	60	REAL(wp), PUBLIC, DIMENSION(jpk) :: & !:
	61	xnumm !: maximum number of particles in aggregates
	62
	63	#endif
	64
	65	#if ! defined key_kriest
	66	REAL(wp), DIMENSION(jpi,jpj,jpk) :: & !:
	67	sinkfer2 !: Big Fe sinking flux
	68	#endif
	69
	70	!!* Substitution
[1808]	71	# include "top_substitute.h90"
[935]	72	!!----------------------------------------------------------------------
	73	!! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
[1152]	74	!! $Id$
[935]	75	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
	76	!!----------------------------------------------------------------------
	77
	78	CONTAINS
	79
	80	#if defined key_kriest
	81
	82	SUBROUTINE p4z_sink ( kt, jnt )
	83	!!---------------------------------------------------------------------
	84	!! * ROUTINE p4z_sink *
	85	!!
	86	!! ** Purpose : Compute vertical flux of particulate matter due to
	87	!! gravitational sinking - Kriest parameterization
	88	!!
	89	!! ** Method : - ???
	90	!!---------------------------------------------------------------------
	91
	92	INTEGER, INTENT(in) :: kt, jnt
	93	INTEGER :: ji, jj, jk
	94	REAL(wp) :: zagg1, zagg2, zagg3, zagg4, zagg5, zaggsi, zaggsh
	95	REAL(wp) :: zagg , zaggdoc, znumdoc
	96	REAL(wp) :: znum , zeps, zfm, zgm, zsm
	97	REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5
	98	REAL(wp) :: zval1, zval2, zval3, zval4
[1457]	99	#if defined key_trc_diaadd
	100	REAL(wp) :: zrfact2
	101	INTEGER :: iksed1
	102	#if defined key_iomput
	103	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d
[935]	104	#endif
[1457]	105	#endif
[935]	106	REAL(wp), DIMENSION(jpi,jpj,jpk) :: znum3d
	107	CHARACTER (len=25) :: charout
	108
	109	!!---------------------------------------------------------------------
	110
[1264]	111	IF( ( kt * jnt ) == nittrc000 ) THEN
	112	CALL p4z_sink_init ! Initialization (first time-step only)
[1736]	113	xstep = rfact2 / rday ! Time step duration for biology
[1264]	114	xstep2 = rfact2 / 2.
	115	ENDIF
[935]	116
	117	! Initialisation of variables used to compute Sinking Speed
	118	! ---------------------------------------------------------
	119
	120	znum3d(:,:,:) = 0.e0
	121	zval1 = 1. + xkr_zeta
	122	zval2 = 1. + xkr_zeta + xkr_eta
	123	zval3 = 1. + xkr_eta
	124
	125	! Computation of the vertical sinking speed : Kriest et Evans, 2000
	126	! -----------------------------------------------------------------
	127
	128	DO jk = 1, jpkm1
	129	DO jj = 1, jpj
	130	DO ji = 1, jpi
	131	IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
	132	znum = trn(ji,jj,jk,jppoc) / ( trn(ji,jj,jk,jpnum) + rtrn ) / xkr_massp
	133	! -------------- To avoid sinking speed over 50 m/day -------
	134	znum = MIN( xnumm(jk), znum )
	135	znum = MAX( 1.1 , znum )
	136	znum3d(ji,jj,jk) = znum
	137	!------------------------------------------------------------
	138	zeps = ( zval1 * znum - 1. )/ ( znum - 1. )
	139	zfm = xkr_frac**( 1. - zeps )
	140	zgm = xkr_frac**( zval1 - zeps )
	141	zdiv = MAX( 1.e-4, ABS( zeps - zval2 ) ) * SIGN( 1., ( zeps - zval2 ) )
	142	zdiv1 = zeps - zval3
	143	wsbio3(ji,jj,jk) = xkr_wsbio_min * ( zeps - zval1 ) / zdiv &
	144	& - xkr_wsbio_max * zgm * xkr_eta / zdiv
	145	wsbio4(ji,jj,jk) = xkr_wsbio_min * ( zeps-1. ) / zdiv1 &
	146	& - xkr_wsbio_max * zfm * xkr_eta / zdiv1
	147	IF( znum == 1.1) wsbio3(ji,jj,jk) = wsbio4(ji,jj,jk)
	148	ENDIF
	149	END DO
	150	END DO
	151	END DO
	152
	153	wscal(:,:,:) = MAX( wsbio3(:,:,:), 50. )
	154
	155
	156	! INITIALIZE TO ZERO ALL THE SINKING ARRAYS
	157	! -----------------------------------------
	158
	159	sinking (:,:,:) = 0.e0
	160	sinking2(:,:,:) = 0.e0
	161	sinkcal (:,:,:) = 0.e0
	162	sinkfer (:,:,:) = 0.e0
	163	sinksil (:,:,:) = 0.e0
	164
	165	! Compute the sedimentation term using p4zsink2 for all
	166	! the sinking particles
	167	! -----------------------------------------------------
	168
	169	CALL p4z_sink2( wsbio3, sinking , jppoc )
	170	CALL p4z_sink2( wsbio4, sinking2, jpnum )
	171	CALL p4z_sink2( wsbio3, sinkfer , jpsfe )
	172	CALL p4z_sink2( wscal , sinksil , jpdsi )
	173	CALL p4z_sink2( wscal , sinkcal , jpcal )
	174
	175	! Exchange between organic matter compartments due to
	176	! coagulation/disaggregation
	177	! ---------------------------------------------------
	178
	179	zval1 = 1. + xkr_zeta
	180	zval2 = 1. + xkr_eta
	181	zval3 = 3. + xkr_eta
	182	zval4 = 4. + xkr_eta
	183
	184	DO jk = 1,jpkm1
	185	DO jj = 1,jpj
	186	DO ji = 1,jpi
	187	IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
	188
	189	znum = trn(ji,jj,jk,jppoc)/(trn(ji,jj,jk,jpnum)+rtrn) / xkr_massp
	190	! -------------- To avoid sinking speed over 50 m/day -------
	191	znum = min(xnumm(jk),znum)
	192	znum = MAX( 1.1,znum)
	193	!------------------------------------------------------------
	194	zeps = ( zval1 * znum - 1.) / ( znum - 1.)
	195	zdiv = MAX( 1.e-4, ABS( zeps - zval3) ) * SIGN( 1., zeps - zval3 )
	196	zdiv1 = MAX( 1.e-4, ABS( zeps - 4. ) ) * SIGN( 1., zeps - 4. )
	197	zdiv2 = zeps - 2.
	198	zdiv3 = zeps - 3.
	199	zdiv4 = zeps - zval2
	200	zdiv5 = 2.* zeps - zval4
	201	zfm = xkr_frac**( 1.- zeps )
	202	zsm = xkr_frac**xkr_eta
	203
	204	! Part I : Coagulation dependant on turbulence
	205	! ----------------------------------------------
	206
	207	zagg1 = ( 0.163 * trn(ji,jj,jk,jpnum)**2 &
	208	& * 2.( (zfm-1.)(zfmxkr_mass_max3-xkr_mass_min*3) &
	209	& * (zeps-1)/zdiv1 + 3.(zfmxkr_mass_max-xkr_mass_min) &
	210	& * (zfmxkr_mass_max2-xkr_mass_min*2) &
	211	& * (zeps-1.)*2/(zdiv2zdiv3)) &
	212	# if defined key_off_degrad
	213	& *facvol(ji,jj,jk) &
	214	# endif
	215	& )
	216
	217	zagg2 = ( 20.163trn(ji,jj,jk,jpnum)*2zfm* &
	218	& ((xkr_mass_max*3+3.(xkr_mass_max**2 &
	219	& xkr_mass_min(zeps-1.)/zdiv2 &
	220	& +xkr_mass_maxxkr_mass_min2(zeps-1.)/zdiv3) &
	221	& +xkr_mass_min*3(zeps-1)/zdiv1) &
	222	& -zfmxkr_mass_max3(1.+3.*((zeps-1.)/ &
	223	& (zeps-2.)+(zeps-1.)/zdiv3)+(zeps-1.)/zdiv1)) &
	224	# if defined key_off_degrad
	225	& *facvol(ji,jj,jk) &
	226	# endif
	227	& )
	228
	229	zagg3 = ( 0.163trn(ji,jj,jk,jpnum)2zfm*28. * xkr_mass_max**3 &
	230	# if defined key_off_degrad
	231	& *facvol(ji,jj,jk) &
	232	# endif
	233	& )
	234
	235	zaggsh = ( zagg1 + zagg2 + zagg3 ) * rfact2 * xdiss(ji,jj,jk) / 1000.
	236
	237	! Aggregation of small into large particles
	238	! Part II : Differential settling
	239	! ----------------------------------------------
	240
	241	zagg4 = ( 2.3.1410.125trn(ji,jj,jk,jpnum)2 &
	242	& xkr_wsbio_min(zeps-1.)*2 &
	243	& (xkr_mass_min2((1.-zsmzfm)/(zdiv3zdiv4) &
	244	& -(1.-zfm)/(zdiv*(zeps-1.)))- &
	245	& ((zfmzfmxkr_mass_max*2zsm-xkr_mass_min**2) &
	246	& xkr_eta)/(zdivzdiv3*zdiv5) ) &
	247	# if defined key_off_degrad
	248	& *facvol(ji,jj,jk) &
	249	# endif
	250	& )
	251
	252	zagg5 = ( 2.3.1410.125trn(ji,jj,jk,jpnum)*2 &
	253	& (zeps-1.)zfm*xkr_wsbio_min &
	254	& (zsm(xkr_mass_min*2-zfmxkr_mass_max**2) &
	255	& /zdiv3-(xkr_mass_min*2-zfmzsmxkr_mass_max*2) &
	256	& /zdiv) &
	257	# if defined key_off_degrad
	258	& *facvol(ji,jj,jk) &
	259	# endif
	260	& )
	261
[1264]	262	zaggsi = ( zagg4 + zagg5 ) * xstep / 10.
[935]	263
	264	zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi )
	265
	266	! Aggregation of DOC to small particles
	267	! --------------------------------------
	268
	269	zaggdoc = ( 0.4 * trn(ji,jj,jk,jpdoc) &
[1264]	270	& + 1018. * trn(ji,jj,jk,jppoc) ) * xstep &
[935]	271	# if defined key_off_degrad
	272	& * facvol(ji,jj,jk) &
	273	# endif
	274	& * xdiss(ji,jj,jk) * trn(ji,jj,jk,jpdoc)
	275
	276	znumdoc = trn(ji,jj,jk,jpnum) / ( trn(ji,jj,jk,jppoc) + rtrn )
	277	tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zaggdoc
	278	tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zaggdoc * znumdoc - zagg
	279	tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc
	280
	281	ENDIF
	282	END DO
	283	END DO
	284	END DO
	285
[1329]	286	#if defined key_trc_diaadd
[935]	287	zrfact2 = 1.e3 * rfact2r
[1457]	288	iksed1 = iksed + 1
	289	# if ! defined key_iomput
	290	trc2d(:,: ,jp_pcs0_2d + 4) = sinking (:,:,iksed1) * zrfact2 * tmask(:,:,1)
	291	trc2d(:,: ,jp_pcs0_2d + 5) = sinking2(:,:,iksed1) * zrfact2 * tmask(:,:,1)
	292	trc2d(:,: ,jp_pcs0_2d + 6) = sinkfer (:,:,iksed1) * zrfact2 * tmask(:,:,1)
	293	trc2d(:,: ,jp_pcs0_2d + 7) = sinksil (:,:,iksed1) * zrfact2 * tmask(:,:,1)
	294	trc2d(:,: ,jp_pcs0_2d + 8) = sinkcal (:,:,iksed1) * zrfact2 * tmask(:,:,1)
	295	trc3d(:,:,:,jp_pcs0_3d + 11) = sinking (:,:,:) * zrfact2 * tmask(:,:,:)
	296	trc3d(:,:,:,jp_pcs0_3d + 12) = sinking2(:,:,:) * zrfact2 * tmask(:,:,:)
	297	trc3d(:,:,:,jp_pcs0_3d + 13) = sinksil (:,:,:) * zrfact2 * tmask(:,:,:)
	298	trc3d(:,:,:,jp_pcs0_3d + 14) = sinkcal (:,:,:) * zrfact2 * tmask(:,:,:)
	299	trc3d(:,:,:,jp_pcs0_3d + 15) = znum3d (:,:,:) * tmask(:,:,:)
	300	trc3d(:,:,:,jp_pcs0_3d + 16) = wsbio3 (:,:,:) * tmask(:,:,:)
	301	trc3d(:,:,:,jp_pcs0_3d + 17) = wsbio4 (:,:,:) * tmask(:,:,:)
	302	#else
[1678]	303	zw3d(:,:,:) = sinking (:,:,:) * zrfact2 * tmask(:,:,:)
	304	IF( jnt == nrdttrc ) CALL iom_put( "PMO" , zw3d )
	305	zw3d(:,:,:) = sinking2(:,:,:) * zrfact2 * tmask(:,:,:)
	306	IF( jnt == nrdttrc ) CALL iom_put( "PMO2", zw3d )
	307	zw3d(:,:,:) = sinkfer (:,:,:) * zrfact2 * tmask(:,:,:)
	308	IF( jnt == nrdttrc ) CALL iom_put( "ExpFe1", zw3d )
	309	zw3d(:,:,:) = sinksil (:,:,:) * zrfact2 * tmask(:,:,:)
	310	IF( jnt == nrdttrc ) CALL iom_put( "ExpSi", zw3d )
	311	zw3d(:,:,:) = sinkcal (:,:,:) * zrfact2 * tmask(:,:,:)
	312	IF( jnt == nrdttrc ) CALL iom_put( "ExpCaCO3", zw3d )
	313	zw3d(:,:,:) = sinking (:,:,:) * zrfact2 * tmask(:,:,:)
[1457]	314	IF( jnt == nrdttrc ) CALL iom_put( "POCFlx", zw3d )
[1678]	315	zw3d(:,:,:) = sinking2(:,:,:) * zrfact2 * tmask(:,:,:)
[1457]	316	IF( jnt == nrdttrc ) CALL iom_put( "GOCFlx", zw3d )
[1678]	317	zw3d(:,:,:) = sinksil (:,:,:) * zrfact2 * tmask(:,:,:)
[1457]	318	IF( jnt == nrdttrc ) CALL iom_put( "SiFlx", zw3d )
[1678]	319	zw3d(:,:,:) = sinkcal (:,:,:) * zrfact2 * tmask(:,:,:)
[1457]	320	IF( jnt == nrdttrc ) CALL iom_put( "CaCO3Flx", zw3d )
[1678]	321	zw3d(:,:,:) = znum3d (:,:,:) * tmask(:,:,:)
[1457]	322	IF( jnt == nrdttrc ) CALL iom_put( "xnum", zw3d )
[1678]	323	zw3d(:,:,:) = wsbio3 (:,:,:) * tmask(:,:,:)
[1457]	324	IF( jnt == nrdttrc ) CALL iom_put( "W1", zw3d )
[1678]	325	zw3d(:,:,:) = wsbio4 (:,:,:) * tmask(:,:,:)
[1457]	326	IF( jnt == nrdttrc ) CALL iom_put( "W2", zw3d )
[1329]	327	# endif
	328
	329	#endif
[935]	330	!
	331	IF(ln_ctl) THEN ! print mean trends (used for debugging)
	332	WRITE(charout, FMT="('sink')")
	333	CALL prt_ctl_trc_info(charout)
	334	CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
	335	ENDIF
	336
	337	END SUBROUTINE p4z_sink
	338
	339	SUBROUTINE p4z_sink_init
	340	!!----------------------------------------------------------------------
	341	!! * ROUTINE p4z_sink_init *
	342	!!
	343	!! ** Purpose : Initialization of sinking parameters
	344	!! Kriest parameterization only
	345	!!
[1119]	346	!! ** Method : Read the nampiskrs namelist and check the parameters
[935]	347	!! called at the first timestep (nittrc000)
	348	!!
[1119]	349	!! ** input : Namelist nampiskrs
[935]	350	!!
	351	!!----------------------------------------------------------------------
	352	INTEGER :: jk, jn, kiter
	353	REAL(wp) :: znum, zdiv
	354	REAL(wp) :: zws, zwr, zwl,wmax, znummax
	355	REAL(wp) :: zmin, zmax, zl, zr, xacc
	356
[1119]	357	NAMELIST/nampiskrs/ xkr_sfact, xkr_stick , &
[935]	358	& xkr_nnano, xkr_ndiat, xkr_nmeso, xkr_naggr
	359
	360	!!----------------------------------------------------------------------
[1119]	361	REWIND( numnat ) ! read nampiskrs
	362	READ ( numnat, nampiskrs )
[935]	363
	364	IF(lwp) THEN
	365	WRITE(numout,*)
[1119]	366	WRITE(numout,*) ' Namelist : nampiskrs'
[935]	367	WRITE(numout,*) ' Sinking factor xkr_sfact = ', xkr_sfact
	368	WRITE(numout,*) ' Stickiness xkr_stick = ', xkr_stick
	369	WRITE(numout,*) ' Nbr of cell in nano size class xkr_nnano = ', xkr_nnano
	370	WRITE(numout,*) ' Nbr of cell in diatoms size class xkr_ndiat = ', xkr_ndiat
	371	WRITE(numout,*) ' Nbr of cell in mesozoo size class xkr_nmeso = ', xkr_nmeso
	372	WRITE(numout,*) ' Nbr of cell in aggregates size class xkr_naggr = ', xkr_naggr
	373	ENDIF
	374
	375
	376	! max and min vertical particle speed
	377	xkr_wsbio_min = xkr_sfact * xkr_mass_min**xkr_eta
	378	xkr_wsbio_max = xkr_sfact * xkr_mass_max**xkr_eta
	379	WRITE(numout,*) ' max and min vertical particle speed ', xkr_wsbio_min, xkr_wsbio_max
	380
	381	!
	382	! effect of the sizes of the different living pools on particle numbers
	383	! nano = 2um-20um -> mean size=6.32 um -> ws=2.596 -> xnum=xnnano=2.337
	384	! diat and microzoo = 10um-200um -> 44.7 -> 8.732 -> xnum=xndiat=3.718
	385	! mesozoo = 200um-2mm -> 632.45 -> 45.14 -> xnum=xnmeso=7.147
	386	! aggregates = 200um-10mm -> 1414 -> 74.34 -> xnum=xnaggr=9.877
	387	! doc aggregates = 1um
	388	! ----------------------------------------------------------
	389
	390	xkr_dnano = 1. / ( xkr_massp * xkr_nnano )
	391	xkr_ddiat = 1. / ( xkr_massp * xkr_ndiat )
	392	xkr_dmeso = 1. / ( xkr_massp * xkr_nmeso )
	393	xkr_daggr = 1. / ( xkr_massp * xkr_naggr )
	394
	395	!!---------------------------------------------------------------------
	396	!! 'key_kriest' ???
	397	!!---------------------------------------------------------------------
	398	! COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED
	399	! Search of the maximum number of particles in aggregates for each k-level.
	400	! Bissection Method
	401	!--------------------------------------------------------------------
	402	WRITE(numout,*)
	403	WRITE(numout,*)' kriest : Compute maximum number of particles in aggregates'
	404
	405	xacc = 0.001
	406	kiter = 50
	407	zmin = 1.10
	408	zmax = xkr_mass_max / xkr_mass_min
	409	xkr_frac = zmax
	410
	411	DO jk = 1,jpk
	412	zl = zmin
	413	zr = zmax
[1736]	414	wmax = 0.5 * fse3t(1,1,jk) * rday / rfact2
[935]	415	zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
	416	znum = zl - 1.
	417	zwl = xkr_wsbio_min * xkr_zeta / zdiv &
	418	& - ( xkr_wsbio_max * xkr_eta * znum * &
	419	& xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
	420	& - wmax
	421
	422	zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
	423	znum = zr - 1.
	424	zwr = xkr_wsbio_min * xkr_zeta / zdiv &
	425	& - ( xkr_wsbio_max * xkr_eta * znum * &
	426	& xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
	427	& - wmax
	428	iflag: DO jn = 1, kiter
	429	IF( zwl == 0.e0 ) THEN
	430	znummax = zl
	431	ELSE IF ( zwr == 0.e0 ) THEN
	432	znummax = zr
	433	ELSE
	434	znummax = ( zr + zl ) / 2.
	435	zdiv = xkr_zeta + xkr_eta - xkr_eta * znummax
	436	znum = znummax - 1.
	437	zws = xkr_wsbio_min * xkr_zeta / zdiv &
	438	& - ( xkr_wsbio_max * xkr_eta * znum * &
	439	& xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
	440	& - wmax
	441	IF( zws * zwl < 0. ) THEN
	442	zr = znummax
	443	ELSE
	444	zl = znummax
	445	ENDIF
	446	zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
	447	znum = zl - 1.
	448	zwl = xkr_wsbio_min * xkr_zeta / zdiv &
	449	& - ( xkr_wsbio_max * xkr_eta * znum * &
	450	& xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
	451	& - wmax
	452
	453	zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
	454	znum = zr - 1.
	455	zwr = xkr_wsbio_min * xkr_zeta / zdiv &
	456	& - ( xkr_wsbio_max * xkr_eta * znum * &
	457	& xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
	458	& - wmax
	459
	460	IF ( ABS ( zws ) <= xacc ) EXIT iflag
	461
	462	ENDIF
	463
	464	END DO iflag
	465
	466	xnumm(jk) = znummax
	467	WRITE(numout,*) ' jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)
	468
	469	END DO
	470
	471	END SUBROUTINE p4z_sink_init
	472
	473	#else
	474
	475	SUBROUTINE p4z_sink ( kt, jnt )
	476	!!---------------------------------------------------------------------
	477	!! * ROUTINE p4z_sink *
	478	!!
	479	!! ** Purpose : Compute vertical flux of particulate matter due to
	480	!! gravitational sinking
	481	!!
	482	!! ** Method : - ???
	483	!!---------------------------------------------------------------------
	484	INTEGER, INTENT(in) :: kt, jnt
	485	INTEGER :: ji, jj, jk
	486	REAL(wp) :: zagg1, zagg2, zagg3, zagg4
	487	REAL(wp) :: zagg , zaggfe, zaggdoc, zaggdoc2
[1264]	488	REAL(wp) :: zfact, zwsmax
[935]	489	#if defined key_trc_dia3d
	490	REAL(wp) :: zrfact2
[1457]	491	INTEGER :: iksed1
[935]	492	#endif
[1457]	493	#if defined key_iomput
[1678]	494	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d
[1457]	495	#endif
[935]	496	CHARACTER (len=25) :: charout
	497	!!---------------------------------------------------------------------
	498
[1264]	499	IF( ( kt * jnt ) == nittrc000 ) THEN
[1736]	500	xstep = rfact2 / rday ! Timestep duration for biology
[1264]	501	xstep2 = rfact2 / 2.
	502	ENDIF
[935]	503
	504	! Sinking speeds of detritus is increased with depth as shown
	505	! by data and from the coagulation theory
	506	! -----------------------------------------------------------
	507	DO jk = 1, jpkm1
	508	DO jj = 1, jpj
	509	DO ji=1,jpi
	510	zfact = MAX( 0., fsdepw(ji,jj,jk+1)-hmld(ji,jj) ) / 4000.
	511	wsbio4(ji,jj,jk) = wsbio2 + ( 200.- wsbio2 ) * zfact
	512	END DO
	513	END DO
	514	END DO
	515
	516	! LIMIT THE VALUES OF THE SINKING SPEEDS
	517	! TO AVOID NUMERICAL INSTABILITIES
	518
	519	wsbio3(:,:,:) = wsbio
	520	!
	521	! OA Below, this is garbage. the ideal would be to find a time-splitting
	522	! OA algorithm that does not increase the computing cost by too much
	523	! OA In ROMS, I have included a time-splitting procedure. But it is
	524	! OA too expensive as the loop is computed globally. Thus, a small e3t
	525	! OA at one place determines the number of subtimesteps globally
	526	! OA AWFULLY EXPENSIVE !! Not able to find a better approach. Damned !!
	527
	528	DO jk = 1,jpkm1
	529	DO jj = 1, jpj
	530	DO ji = 1, jpi
[1264]	531	zwsmax = 0.8 * fse3t(ji,jj,jk) / xstep
[935]	532	wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax )
	533	wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax )
	534	END DO
	535	END DO
	536	END DO
	537
	538	wscal(:,:,:) = wsbio4(:,:,:)
	539
	540	! INITIALIZE TO ZERO ALL THE SINKING ARRAYS
	541	! -----------------------------------------
	542
	543	sinking (:,:,:) = 0.e0
	544	sinking2(:,:,:) = 0.e0
	545	sinkcal (:,:,:) = 0.e0
	546	sinkfer (:,:,:) = 0.e0
	547	sinksil (:,:,:) = 0.e0
	548	sinkfer2(:,:,:) = 0.e0
	549
	550	! Compute the sedimentation term using p4zsink2 for all
	551	! the sinking particles
	552	! -----------------------------------------------------
	553
	554	CALL p4z_sink2( wsbio3, sinking , jppoc )
	555	CALL p4z_sink2( wsbio3, sinkfer , jpsfe )
	556	CALL p4z_sink2( wsbio4, sinking2, jpgoc )
	557	CALL p4z_sink2( wsbio4, sinkfer2, jpbfe )
	558	CALL p4z_sink2( wsbio4, sinksil , jpdsi )
	559	CALL p4z_sink2( wscal , sinkcal , jpcal )
	560
	561	! Exchange between organic matter compartments due to
	562	! coagulation/disaggregation
	563	! ---------------------------------------------------
	564
	565	DO jk = 1, jpkm1
	566	DO jj = 1, jpj
	567	DO ji = 1, jpi
[1264]	568	zfact = xstep * xdiss(ji,jj,jk)
[1457]	569	! Part I : Coagulation dependent on turbulence
[935]	570	# if defined key_off_degrad
	571	zagg1 = 940.* zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc) * facvol(ji,jj,jk)
[1457]	572	zagg2 = 1.054e4 * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc) * facvol(ji,jj,jk)
[935]	573	# else
	574	zagg1 = 940.* zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
	575	zagg2 = 1.054e4 * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
	576	# endif
	577
[1457]	578	! Part II : Differential settling
[935]	579
[1457]	580	! Aggregation of small into large particles
[935]	581	# if defined key_off_degrad
[1264]	582	zagg3 = 0.66 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc) * facvol(ji,jj,jk)
[1457]	583	zagg4 = 0.e0 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc) * facvol(ji,jj,jk)
[935]	584	# else
[1264]	585	zagg3 = 0.66 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
	586	zagg4 = 0.e0 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
[935]	587	# endif
	588	zagg = zagg1 + zagg2 + zagg3 + zagg4
	589	zaggfe = zagg * trn(ji,jj,jk,jpsfe) / ( trn(ji,jj,jk,jppoc) + rtrn )
	590
[1457]	591	! Aggregation of DOC to small particles
	592	#if defined key_off_degrad
[935]	593	zaggdoc = ( 80.* trn(ji,jj,jk,jpdoc) + 698. * trn(ji,jj,jk,jppoc) ) &
[1457]	594	& * facvol(ji,jj,jk) * zfact * trn(ji,jj,jk,jpdoc)
[935]	595	zaggdoc2 = 1.05e4 * zfact * trn(ji,jj,jk,jpgoc) &
[1457]	596	& * facvol(ji,jj,jk) * trn(ji,jj,jk,jpdoc)
	597	#else
	598	zaggdoc = ( 80.* trn(ji,jj,jk,jpdoc) + 698. * trn(ji,jj,jk,jppoc) ) &
	599	& * zfact * trn(ji,jj,jk,jpdoc)
	600	zaggdoc2 = 1.05e4 * zfact * trn(ji,jj,jk,jpgoc) * trn(ji,jj,jk,jpdoc)
	601	#endif
	602	! Update the trends
[935]	603	tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zagg + zaggdoc
	604	tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zagg + zaggdoc2
	605	tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zaggfe
	606	tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggfe
	607	tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2
[1457]	608	!
[935]	609	END DO
	610	END DO
	611	END DO
	612
[1457]	613	#if defined key_trc_diaadd
[935]	614	zrfact2 = 1.e3 * rfact2r
[1457]	615	iksed1 = iksed + 1
	616	# if ! defined key_iomput
	617	trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,iksed1) * zrfact2 * tmask(:,:,1)
	618	trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,iksed1) * zrfact2 * tmask(:,:,1)
	619	trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,iksed1) * zrfact2 * tmask(:,:,1)
	620	trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,iksed1) * zrfact2 * tmask(:,:,1)
	621	trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,iksed1) * zrfact2 * tmask(:,:,1)
	622	trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,iksed1) * zrfact2 * tmask(:,:,1)
	623	# else
[1678]	624	zw3d(:,:,:) = sinking (:,:,:) * zrfact2 * tmask(:,:,:)
	625	IF( jnt == nrdttrc ) CALL iom_put( "ExpPOC" , zw3d )
	626	zw3d(:,:,:) = sinking2(:,:,:) * zrfact2 * tmask(:,:,:)
	627	IF( jnt == nrdttrc ) CALL iom_put( "ExpGOC", zw3d )
	628	zw3d(:,:,:) = sinkfer (:,:,:) * zrfact2 * tmask(:,:,:)
	629	IF( jnt == nrdttrc ) CALL iom_put( "ExpFe1", zw3d )
	630	zw3d(:,:,:) = sinkfer2(:,:,:) * zrfact2 * tmask(:,:,:)
	631	IF( jnt == nrdttrc ) CALL iom_put( "ExpFe2", zw3d )
	632	zw3d(:,:,:) = sinksil (:,:,:) * zrfact2 * tmask(:,:,:)
	633	IF( jnt == nrdttrc ) CALL iom_put( "ExpSi", zw3d )
	634	zw3d(:,:,:) = sinkcal (:,:,:) * zrfact2 * tmask(:,:,:)
	635	IF( jnt == nrdttrc ) CALL iom_put( "Expcal", zw3d )
[1457]	636	# endif
	637	#endif
[935]	638	!
	639	IF(ln_ctl) THEN ! print mean trends (used for debugging)
	640	WRITE(charout, FMT="('sink')")
	641	CALL prt_ctl_trc_info(charout)
	642	CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
	643	ENDIF
	644
	645	END SUBROUTINE p4z_sink
	646
	647	#endif
	648
[1073]	649	SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra )
	650	!!---------------------------------------------------------------------
	651	!! * ROUTINE p4z_sink2 *
	652	!!
	653	!! ** Purpose : Compute the sedimentation terms for the various sinking
	654	!! particles. The scheme used to compute the trends is based
	655	!! on MUSCL.
	656	!!
	657	!! ** Method : - this ROUTINE compute not exactly the advection but the
	658	!! transport term, i.e. div(u*tra).
	659	!!---------------------------------------------------------------------
	660	INTEGER , INTENT(in ) :: jp_tra ! tracer index index
	661	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwsink ! sinking speed
	662	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: psinkflx ! sinking fluxe
	663	!!
	664	INTEGER :: ji, jj, jk, jn
[1264]	665	REAL(wp) :: zigma,zew,zign, zflx
[1073]	666	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztraz, zakz
	667	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwsink2
	668	!!---------------------------------------------------------------------
	669
	670
	671	ztraz(:,:,:) = 0.e0
	672	zakz (:,:,:) = 0.e0
	673
	674	DO jk = 1, jpkm1
	675	# if defined key_off_degrad
[1736]	676	zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1) * facvol(:,:,jk)
[1073]	677	# else
[1736]	678	zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1)
[1073]	679	# endif
	680	END DO
	681	zwsink2(:,:,1) = 0.e0
	682
	683
	684	! Vertical advective flux
	685	DO jn = 1, 2
	686	! first guess of the slopes interior values
	687	DO jk = 2, jpkm1
	688	ztraz(:,:,jk) = ( trn(:,:,jk-1,jp_tra) - trn(:,:,jk,jp_tra) ) * tmask(:,:,jk)
	689	END DO
	690	ztraz(:,:,1 ) = 0.0
	691	ztraz(:,:,jpk) = 0.0
	692
	693	! slopes
	694	DO jk = 2, jpkm1
	695	DO jj = 1,jpj
	696	DO ji = 1, jpi
	697	zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
	698	zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
	699	END DO
	700	END DO
	701	END DO
	702
	703	! Slopes limitation
	704	DO jk = 2, jpkm1
	705	DO jj = 1, jpj
	706	DO ji = 1, jpi
	707	zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) * &
	708	& MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
	709	END DO
	710	END DO
	711	END DO
	712
	713	! vertical advective flux
	714	DO jk = 1, jpkm1
	715	DO jj = 1, jpj
	716	DO ji = 1, jpi
[1264]	717	zigma = zwsink2(ji,jj,jk+1) * xstep2 / fse3w(ji,jj,jk+1)
[1073]	718	zew = zwsink2(ji,jj,jk+1)
[1264]	719	psinkflx(ji,jj,jk+1) = -zew * ( trn(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * xstep2
[1073]	720	END DO
	721	END DO
	722	END DO
	723	!
	724	! Boundary conditions
	725	psinkflx(:,:,1 ) = 0.e0
	726	psinkflx(:,:,jpk) = 0.e0
	727
	728	DO jk=1,jpkm1
	729	DO jj = 1,jpj
	730	DO ji = 1, jpi
	731	zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
	732	trn(ji,jj,jk,jp_tra) = trn(ji,jj,jk,jp_tra) + zflx
	733	END DO
	734	END DO
	735	END DO
	736
	737	ENDDO
	738
	739	DO jk=1,jpkm1
	740	DO jj = 1,jpj
	741	DO ji = 1, jpi
	742	zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
	743	trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + 2. * zflx
	744	END DO
	745	END DO
	746	END DO
	747
	748	trn(:,:,:,jp_tra) = trb(:,:,:,jp_tra)
	749	psinkflx(:,:,:) = 2. * psinkflx(:,:,:)
	750
	751	!
	752	END SUBROUTINE p4z_sink2
	753
[935]	754	#else
	755	!!======================================================================
	756	!! Dummy module : No PISCES bio-model
	757	!!======================================================================
	758	CONTAINS
	759	SUBROUTINE p4z_sink ! Empty routine
	760	END SUBROUTINE p4z_sink
	761	#endif
	762
	763	!!======================================================================
	764	END MODULE p4zsink

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

Download in other formats: