!!---------------------------------------------------------------------- !! *** trcini.pisces.h90 *** !!---------------------------------------------------------------------- # include "domzgr_substitute.h90" # include "passivetrc_substitute.h90" CONTAINS SUBROUTINE trc_ini !!----------------------------------------------------------------- !! !! *** ROUTINE trc_ini *** !! !! !! Purpose : !! --------- !! Initialisation of PISCES biological and chemical variables !! !! INPUT : !! ----- !! common !! all the common defined in opa !! !! !! OUTPUT : : no !! ------ !! !! EXTERNAL : !! ---------- !! p4zche !! !! MODIFICATIONS: !! -------------- !! original : 1988-07 E. MAIER-REIMER MPI HAMBURG !! additions : 1999-10 O. Aumont and C. Le Quere !! additions : 2002 O. Aumont (PISCES) !! 03-2005 O. Aumont and A. El Moussaoui F90 !!---------------------------------------------------------------------- !! TOP 1.0 , LOCEAN-IPSL (2005) !! $Header$ !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt !!---------------------------------------------------------------------- !!Module used USE iom !! local declarations !! ================== INTEGER :: ji,jj,jk INTEGER :: ichl,iband,jm INTEGER , PARAMETER :: jpmois = 12, jpan = 1 REAL(wp) :: ztoto,expide,denitide,ztra,zmaskt REAL(wp) , DIMENSION (jpi,jpj) :: riverdoc,river,ndepo REAL(wp) , DIMENSION (jpi,jpj,jpk) :: cmask INTEGER :: numriv,numdust,numbath,numdep INTEGER :: numlight #if defined key_trc_kriest REAL(wp) :: & znum, zdiv, & zws,zwr, zwl,wmax, xnummax, & zmin, zmax, zl, zr, xacc INTEGER :: jn, kiter #endif !! 1. initialization !! ----------------- !! computation of the record length for direct access FILE !! this length depend of 512 for the t3d machine !! rfact = rdttra(1) * float(ndttrc) rfactr = 1./rfact IF(lwp) WRITE(numout,*) ' Tracer time step=',rfact,' rdt=',rdt rfact2= rfact / float(nrdttrc) rfact2r = 1./rfact2 IF(lwp) write(numout,*) ' Biology time step=',rfact2 !! INITIALISE DUST INPUT FROM ATMOSPHERE !! ------------------------------------- IF ( bdustfer ) THEN IF(lwp) WRITE(numout,*) ' Initialize dust input from atmosphere ' CALL iom_open ( 'dust.orca.nc', numdust ) DO jm = 1, jpmois CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,jm), jm ) ENDDO CALL iom_close( numdust ) ELSE dustmo(:,:,:) = 0. ENDIF !! INITIALISE THE NUTRIENT INPUT BY RIVERS !! --------------------------------------- IF ( briver ) THEN IF(lwp) WRITE(numout,*) ' Initialize the nutrient input by rivers ' CALL iom_open ( 'river.orca.nc', numriv ) CALL iom_get ( numriv, jpdom_data, 'riverdic', river (:,:), jpan ) CALL iom_get ( numriv, jpdom_data, 'riverdoc', riverdoc(:,:), jpan ) CALL iom_close( numriv ) ELSE river (:,:) = 0. riverdoc(:,:) = 0. endif !! INITIALISE THE N INPUT BY DUST !! --------------------------------------- IF ( bndepo ) THEN IF(lwp) WRITE(numout,*) ' Initialize the nutrient input by dust ' CALL iom_open ( 'ndeposition.orca.nc', numdep ) CALL iom_get ( numdep, jpdom_data, 'ndep', ndepo(:,:), jpan ) CALL iom_close( numdep ) ELSE ndepo(:,:) = 0. ENDIF !! Computation of the coastal mask. !! Computation of an island mask to enhance coastal supply of iron !! --------------------------------------------------------------- IF ( bsedinput ) THEN IF(lwp) WRITE(numout,*) ' Computation of an island mask to enhance coastal supply of iron ' CALL iom_open ( 'bathy.orca.nc', numbath ) CALL iom_get ( numbath, jpdom_data, 'bathy', cmask(:,:,:), jpan ) DO jk = 1, 5 DO jj = 2, jpjm1 DO ji = 2, jpim1 IF ( tmask(ji,jj,jk) /= 0. ) THEN zmaskt = tmask(ji+1,jj,jk) * tmask(ji-1,jj,jk) * tmask(ji,jj+1,jk) & & * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1) IF ( zmaskt == 0. ) THEN cmask(ji,jj,jk ) = 0.1 ENDIF ENDIF END DO END DO END DO DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi expide = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) ) denitide = -0.9543 + 0.7662 * LOG( expide ) - 0.235 * LOG( expide )**2 cmask(ji,jj,jk) = cmask(ji,jj,jk) * MIN( 1., EXP( denitide ) / 0.5 ) END DO END DO END DO CALL iom_close( numbath ) ELSE cmask(:,:,:) = 0. ENDIF ! Lateral boundary conditions on ( avt, en ) (sign unchanged) CALL lbc_lnk( cmask , 'T', 1. ) !! Computation of the total atmospheric supply of Si !! ------------------------------------------------- sumdepsi = 0. DO jm = 1, jpmois DO jj = 2, jpjm1 DO ji = 2, jpim1 sumdepsi = sumdepsi + dustmo(ji,jj,jm)/(12.*rmoss)*8.8 & *0.075/28.1*e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,1) END DO END DO END DO IF( lk_mpp ) CALL mpp_sum( sumdepsi ) ! sum over the global domain !! COMPUTATION OF THE N/P RELEASE DUE TO COASTAL RIVERS !! COMPUTATION OF THE Si RELEASE DUE TO COASTAL RIVERS !! --------------------------------------------------- DO jj=1,jpj DO ji=1,jpi cotdep(ji,jj)=river(ji,jj)*1E9/(12.*raass & *e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,1)+rtrn)*tmask(ji,jj,1) rivinp(ji,jj)=(river(ji,jj)+riverdoc(ji,jj))*1E9 & /(31.6*raass*e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,1)+rtrn) & *tmask(ji,jj,1) nitdep(ji,jj)=7.6*ndepo(ji,jj)*tmask(ji,jj,1)/(14E6*raass & *fse3t(ji,jj,1)+rtrn) END DO END DO ! Lateral boundary conditions on ( cotdep, rivinp, nitdep ) (sign unchanged) CALL lbc_lnk( cotdep , 'T', 1. ) ; CALL lbc_lnk( rivinp , 'T', 1. ) ; CALL lbc_lnk( nitdep , 'T', 1. ) rivpo4input=0. rivalkinput=0. nitdepinput=0. DO jj=2,jpjm1 DO ji=2,jpim1 rivpo4input=rivpo4input+rivinp(ji,jj)*(e1t(ji,jj)*e2t(ji,jj) & *fse3t(ji,jj,1))*tmask(ji,jj,1)*raass rivalkinput=rivalkinput+cotdep(ji,jj)*(e1t(ji,jj)*e2t(ji,jj) & *fse3t(ji,jj,1))*tmask(ji,jj,1)*raass nitdepinput=nitdepinput+nitdep(ji,jj)*(e1t(ji,jj)*e2t(ji,jj) & *fse3t(ji,jj,1))*tmask(ji,jj,1)*raass END DO END DO IF( lk_mpp ) THEN CALL mpp_sum( rivpo4input ) ! sum over the global domain CALL mpp_sum( rivalkinput ) ! sum over the global domain CALL mpp_sum( nitdepinput ) ! sum over the global domain ENDIF !! Coastal supply of iron !! ---------------------- DO jk=1,jpkm1 ironsed(:,:,jk)=sedfeinput*cmask(:,:,jk) & /(fse3t(:,:,jk)*rjjss) END DO ! Lateral boundary conditions on ( ironsed ) (sign unchanged) CALL lbc_lnk( ironsed , 'T', 1. ) #if defined key_trc_kriest !---------------------------------------------------------------------- ! COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED ! Search of the maximum number of particles in aggregates for each k-level. ! Bissection Method !-------------------------------------------------------------------- WRITE(numout,*)' ' WRITE(numout,*)'Compute maximum number of particles in aggregates' WRITE(numout,*)' ' xacc = 0.001 kiter = 50 zmin = 1.10 zmax = xkr_mass_max/xkr_mass_min xkr_frac = zmax DO jk =1,jpk zl = zmin zr = zmax wmax = 0.5 * fse3t(1,1,jk) * rjjss / rfact2 zdiv = xkr_zeta + xkr_eta - xkr_eta * zl znum = zl - 1. zwl = xkr_wsbio_min * xkr_zeta / zdiv & & - ( xkr_wsbio_max * xkr_eta * znum * & & xkr_frac**( -xkr_zeta / znum ) / zdiv ) & & - wmax zdiv = xkr_zeta + xkr_eta - xkr_eta * zr znum = zr - 1. zwr = xkr_wsbio_min * xkr_zeta / zdiv & & - ( xkr_wsbio_max * xkr_eta * znum * & & xkr_frac**( -xkr_zeta / znum ) / zdiv ) & & - wmax iflag: DO jn = 1, kiter IF( zwl == 0. ) THEN xnummax = zl ELSE IF ( zwr == 0. ) THEN xnummax = zr ELSE xnummax = ( zr + zl ) / 2. zdiv = xkr_zeta + xkr_eta - xkr_eta * xnummax znum = xnummax - 1. zws = xkr_wsbio_min * xkr_zeta / zdiv & & - ( xkr_wsbio_max * xkr_eta * znum * & & xkr_frac**( -xkr_zeta / znum ) / zdiv ) & & - wmax IF( zws * zwl < 0. ) THEN zr = xnummax ELSE zl = xnummax ENDIF zdiv = xkr_zeta + xkr_eta - xkr_eta * zl znum = zl - 1. zwl = xkr_wsbio_min * xkr_zeta / zdiv & & - ( xkr_wsbio_max * xkr_eta * znum * & & xkr_frac**( -xkr_zeta / znum ) / zdiv ) & & - wmax zdiv = xkr_zeta + xkr_eta - xkr_eta * zr znum = zr - 1. zwr = xkr_wsbio_min * xkr_zeta / zdiv & & - ( xkr_wsbio_max * xkr_eta * znum * & & xkr_frac**( -xkr_zeta / znum ) / zdiv ) & & - wmax IF ( ABS ( zws ) <= xacc ) EXIT iflag ENDIF ENDDO iflag xnumm(jk) = xnummax WRITE(numout,*) 'jk = ',jk,' wmax = ',wmax,' xnum max = ',xnumm(jk) END DO WRITE(numout,*) '------------------------------------' #endif !!---------------------------------------------------------------------- !! !! Initialize biological variables !! !!---------------------------------------------------------------------- !! Set biological ratios !! --------------------- rno3 = (16.+2.)/122. po4r = 1./122. o2nit = 32./122. rdenit = 97.6/16. o2ut = 140./122. !!---------------------------------------------------------------------- !! !! Initialize chemical variables !! !!---------------------------------------------------------------------- !! set pre-industrial atmospheric [co2] (ppm) and o2/n2 ratio !! ---------------------------------------------------------- atcox = 0.20946 !! Set lower/upper limits for temperature and salinity !! --------------------------------------------------- salchl = 1./1.80655 calcon = 1.03E-2 !! Set coefficients for apparent solubility equilibrium of calcite !! Millero et al. 1995 from Mucci 1983 !! -------------------------------------------------------------- akcc1 = -171.9065 akcc2 = -0.077993 akcc3 = 2839.319 akcc4 = 71.595 akcc5 = -0.77712 akcc6 = 0.0028426 akcc7 = 178.34 akcc8 = -0.07711 akcc9 = 0.0041249 !! Set coefficients for seawater pressure correction !! ------------------------------------------------- devk1(1) = -25.5 devk2(1) = 0.1271 devk3(1) = 0. devk4(1) = -3.08E-3 devk5(1) = 0.0877E-3 !! devk1(2) = -15.82 devk2(2) = -0.0219 devk3(2) = 0. devk4(2) = 1.13E-3 devk5(2) = -0.1475E-3 !! devk1(3) = -29.48 devk2(3) = 0.1622 devk3(3) = 2.608E-3 devk4(3) = -2.84E-3 devk5(3) = 0. !! devk1(4) = -14.51 devk2(4) = 0.1211 devk3(4) = -0.321E-3 devk4(4) = -2.67E-3 devk5(4) = 0.0427E-3 !! devk1(5) = -23.12 devk2(5) = 0.1758 devk3(5) = -2.647E-3 devk4(5) = -5.15E-3 devk5(5) = 0.09E-3 !! devk1(6) = -26.57 devk2(6) = 0.2020 devk3(6) = -3.042E-3 devk4(6) = -4.08E-3 devk5(6) = 0.0714E-3 !! devk1(7) = -25.60 devk2(7) = 0.2324 devk3(7) = -3.6246E-3 devk4(7) = -5.13E-3 devk5(7) = 0.0794E-3 !! !! For calcite with Edmond and Gieske 1970 !! devkst = 0.23 !! devks = 35.4 !! Millero 95 takes this depth dependance for calcite devk1(8) = -48.76 devk2(8) = 0.5304 devk3(8) = 0. devk4(8) = -11.76E-3 devk5(8) = 0.3692E-3 !! !! Coefficients for sulfate and fluoride devk1(9) = -18.03 devk2(9) = 0.0466 devk3(9) = 0.316E-3 devk4(9) = -4.53E-3 devk5(9) = 0.09E-3 devk1(10) = -9.78 devk2(10) = -0.0090 devk3(10) = -0.942E-3 devk4(10) = -3.91E-3 devk5(10) = 0.054E-3 !! Set universal gas constants !! --------------------------- rgas = 83.143 oxyco = 1./22.4144 !! Set boron constants !! ------------------- bor1 = 0.00023 bor2 = 1./10.82 !! Set volumetric solubility constants for co2 in ml/l (Weiss, 1974) !! ----------------------------------------------------------------- c00 = -60.2409 c01 = 93.4517 c02 = 23.3585 c03 = 0.023517 c04 = -0.023656 c05 = 0.0047036 ca0 = -162.8301 ca1 = 218.2968 ca2 = 90.9241 ca3 = -1.47696 ca4 = 0.025695 ca5 = -0.025225 ca6 = 0.0049867 !! Set coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970) !! --------------------------------------------------------------------- c10 = -3670.7 c11 = 62.008 c12 = -9.7944 c13 = 0.0118 c14 = -0.000116 !! Set coeff. for 2. dissoc. of carbonic acid (Edmond and Gieskes, 1970) !! --------------------------------------------------------------------- c20 = -1394.7 c21 = -4.777 c22 = 0.0184 c23 = -0.000118 !! Set constants for calculate concentrations for sulfate and fluoride !! sulfates (Morris & Riley 1966) !!---------------------------------------------------------------------- st1 = 0.14 st2 = 1./96.062 !! fluoride !!------------ ft1 = 0.000067 ft2 = 1./18.9984 !! sulfates (Dickson 1990 change to mol:kg soln, idem OCMIP) !!---------------------------------------------------------- ks0 = 141.328 ks1 = -4276.1 ks2 = -23.093 ks3 = -13856 ks4 = 324.57 ks5 = -47.986 ks6 = 35474 ks7 = -771.54 ks8 = 114.723 ks9 = -2698 ks10 = 1776 ks11 = 1. ks12 = -0.001005 !! fluorides (Dickson & Riley 1979 change to mol/kg soln) !!------------------------------------------------------- kf0 = -12.641 kf1 = 1590.2 kf2 = 1.525 kf3 = 1.0 kf4 = -0.001005 !! !! Set coeff. for 1. dissoc. of boric acid (Edmond and Gieskes, 1970) !! ------------------------------------------------------------------ cb0 = -8966.90 cb1 = -2890.53 cb2 = -77.942 cb3 = 1.728 cb4 = -0.0996 cb5 = 148.0248 cb6 = 137.1942 cb7 = 1.62142 cb8 = -24.4344 cb9 = -25.085 cb10 = -0.2474 cb11 = 0.053105 !! Set coeff. for dissoc. of water (Dickson and Riley, 1979, !! eq. 7, coefficient cw2 corrected from 0.9415 to 0.09415 !! after pers. commun. to B. Bacastow, 1988) !! --------------------------------------------------------- cw0 = -13847.26 cw1 = 148.9652 cw2 = -23.6521 cw3 = 118.67 cw4 = -5.977 cw5 = 1.0495 cw6 = -0.01615 !! Set coeff. for dissoc. of phosphate (Millero (1974) !! --------------------------------------------------- cp10 = 115.54 cp11 = -4576.752 cp12 = -18.453 cp13 = -106.736 cp14 = 0.69171 cp15 = -0.65643 cp16 = -0.01844 cp20 = 172.1033 cp21 = -8814.715 cp22 = -27.927 cp23 = -160.340 cp24 = 1.3566 cp25 = 0.37335 cp26 = -0.05778 cp30 = -18.126 cp31 = -3070.75 cp32 = 17.27039 cp33 = 2.81197 cp34 = -44.99486 cp35 = -0.09984 !! Set coeff. for dissoc. of phosphate (Millero (1974) !! --------------------------------------------------- cs10 = 117.385 cs11 = -8904.2 cs12 = -19.334 cs13 = -458.79 cs14 = 3.5913 cs15 = 188.74 cs16 = -1.5998 cs17 = -12.1652 cs18 = 0.07871 cs19 = 0. cs20 = 1. cs21 = -0.001005 !! Set volumetric solubility constants for o2 in ml/l (Weiss, 1970) !! ---------------------------------------------------------------- ox0 = -58.3877 ox1 = 85.8079 ox2 = 23.8439 ox3 = -0.034892 ox4 = 0.015568 ox5 = -0.0019387 !! FROM THE NEW BIOOPTIC MODEL PROPOSED JM ANDRE, WE READ HERE !! A PRECOMPUTED ARRAY CORRESPONDING TO THE ATTENUATION COEFFICIENT CALL ctlopn( numlight, 'kRGB61.txt', 'OLD', 'FORMATTED', 'SEQUENTIAL', & & 1, numout, .TRUE., 1 ) DO ichl = 1,61 READ(numlight,*) ztoto,(xkrgb(iband,ichl),iband = 1,3) END DO CLOSE(numlight) !! Call p4zche to initialize the chemical constants !! ------------------------------------------------ CALL p4zche !! !! Initialize a counter for the computation of chemistry !! ndayflxtr=0 IF(lwp) WRITE(numout,*) ' Initialisation of PISCES done' END SUBROUTINE trc_ini