#if defined(ROW_LAND) #define SEA_P .true. #define SEA_U .true. #define SEA_V .true. #elif defined(ROW_ALLSEA) #define SEA_P allip(j).or.ip(i,j).ne.0 #define SEA_U alliu(j).or.iu(i,j).ne.0 #define SEA_V alliv(j).or.iv(i,j).ne.0 #else #define SEA_P ip(i,j).ne.0 #define SEA_U iu(i,j).ne.0 #define SEA_V iv(i,j).ne.0 #endif subroutine initrc(mnth) use mod_xc ! HYCOM communication interface use mod_cb_arrays ! HYCOM saved arrays use mod_pipe ! HYCOM debugging interface implicit none ! integer mnth ! ! --- -------------------------- ! --- initializatize all tracers ! --- -------------------------- ! logical lpipe_initrc parameter (lpipe_initrc=.false.) ! character ptxt*12,cformat*99 integer i,ibio,nbio,j,k,ktr real bio_n,bio_p,zk real pwij(kk+1),trwij(kk,ntracr), & prij(kk+1),trcij(kk,ntracr) real chl,swfrac(kdm+1) ! if (ntracr.eq.0) then return ! no tracer endif ! ! --- expand trcflg to allow for number of biology fields. ! nbio = 0 ibio = 0 do ktr= 1,ntracr+1 if (ktr.ne.ntracr+1 .and. & trcflg(min(ktr,ntracr)).eq.9) then if (ibio.eq.0) then !start biology ibio = ktr endif elseif (ibio.ne.0) then !end biology nbio = ktr-ibio if (nbio.eq.3) then ! --- Franks NPZ. trcflg(ibio) = 903 trcflg(ibio+1) = -903 trcflg(ibio+2) = -903 ibio = 0 elseif (nbio.eq.3) then ! --- Two Franks NPZ. trcflg(ibio) = 903 trcflg(ibio+1) = -903 trcflg(ibio+2) = -903 trcflg(ibio+3) = 903 trcflg(ibio+4) = -903 trcflg(ibio+5) = -903 ibio = 0 elseif (nbio.eq.4) then ! --- Lima/Idrisi NPZD. trcflg(ibio) = 904 trcflg(ibio+1) = -904 trcflg(ibio+2) = -904 trcflg(ibio+3) = -904 ibio = 0 elseif (nbio.eq.7) then ! --- Lima/Idrisi NPZD and Franks NPZ. trcflg(ibio) = 904 trcflg(ibio+1) = -904 trcflg(ibio+2) = -904 trcflg(ibio+3) = -904 trcflg(ibio+4) = 903 trcflg(ibio+5) = -903 trcflg(ibio+6) = -903 ibio = 0 elseif (nbio.eq.8) then ! --- Two Lima/Idrisi NPZD. trcflg(ibio) = 904 trcflg(ibio+1) = -904 trcflg(ibio+2) = -904 trcflg(ibio+3) = -904 trcflg(ibio+4) = 904 trcflg(ibio+5) = -904 trcflg(ibio+6) = -904 trcflg(ibio+7) = -904 ibio = 0 elseif (nbio.eq.9) then ! --- Chai 9-component. ! trcflg(ibio) = 909 ! trcflg(ibio+1) = -909 ! trcflg(ibio+2) = -909 ! trcflg(ibio+3) = -909 ! trcflg(ibio+4) = -909 ! trcflg(ibio+5) = -909 ! trcflg(ibio+6) = -909 ! trcflg(ibio+7) = -909 ! trcflg(ibio+8) = -909 ! ibio = 0 ! --- not yet implemented if (mnproc.eq.1) then write(lp,'(/ 3a /)') & 'error - trcflg=9 (standard biology) configured', & ' with 9 consecutive tracers, but Chai scheme is', & ' not yet implemented' call flush(lp) endif !1st tile call xcstop('(trcini)') stop '(trcini)' else ! --- unknown standard biology. if (mnproc.eq.1) then write(lp,'(/ 2a,i3 /)') & 'error - trcflg=9 (standard biology) expects', & ' 3/4/6/7/8 consecutive tracers but have',nbio ! & ' 3/4/6/7/8/9 consecutive tracers but have',nbio call flush(lp) endif !1st tile call xcstop('(trcini)') stop '(trcini)' endif endif enddo ! if (mnproc.eq.1) then write(lp,*) do k= 1,ntracr write(lp,'(a,i3,i6)') 'initrc: k,trcflg =',k,trcflg(k) enddo write(lp,*) endif !1st tile ! if (nbio.gt.0) then ! ! --- input bio-tracer parameters. ! --- note that multiple sets of bio-tracers are allowed, ! --- each is read from tracer.input in tracer order. ! open(unit=uoff+99,file=trim(flnminp)//'tracer.input') do ktr= 1,ntracr if (trcflg(ktr).eq.903) then ! --- NPZ call trcupd_903(1,2, -ktr) elseif (trcflg(ktr).eq.904) then ! --- NPZD call trcupd_904(1,2, -ktr) ! elseif (trcflg(ktr).eq.909) then ! --- Chai 9-component. ! call trcupd_909(1,2, -ktr) endif enddo close(unit=uoff+99) endif ! if (trcrin) then return ! tracer from restart endif ! if (iniflg.eq.2) then ! use climatology call rdrlax(mnth,1) !$OMP PARALLEL DO PRIVATE(j,i,k,ktr,pwij,trwij,prij,trcij) & !$OMP SCHEDULE(STATIC,jblk) do j=1,jj do i=1,ii if (SEA_P) then prij(1)=0.0 do k=1,kk prij(k+1)=prij(k)+dp(i,j,k,1) pwij(k) =pwall(i,j,k,1) do ktr= 1,ntracr trwij(k,ktr)=trwall(i,j,k,1,ktr) enddo !ktr enddo !k pwij(kk+1)=prij(kk+1) ! call plctrc(trwij,pwij,kk,ntracr, ! & trcij,prij,kk ) call plmtrc(trwij,pwij,kk,ntracr, & trcij,prij,kk ) do k=1,kk do ktr= 1,ntracr tracer(i,j,k,1,ktr)=trcij(k,ktr) tracer(i,j,k,2,ktr)=trcij(k,ktr) enddo !ktr enddo !k endif !ip enddo !i enddo !j !$OMP END PARALLEL DO else ! analytic inititalization !$OMP PARALLEL DO PRIVATE(j,i,k,ktr,prij,chl,swfrac,zk,bio_n,bio_p) & !$OMP SCHEDULE(STATIC,jblk) do j=1,jj do i=1,ii if (SEA_P) then prij(1)=0.0 do k=1,kk prij(k+1)=prij(k)+dp(i,j,k,1) enddo !k do ktr= 1,ntracr if (trcflg(ktr).eq.1) then !need the euphotic zone if (jerlv0.le.0) then chl = akpar(i,j,lk0)*wk0+akpar(i,j,lk1)*wk1 & +akpar(i,j,lk2)*wk2+akpar(i,j,lk3)*wk3 endif call swfrac_ij(chl,prij,kk+1,qonem, & !oneta might not be available jerlov(i,j),swfrac) exit !only calculate swfrac once endif !trcflg==1 enddo !ktr do k=1,kk do ktr= 1,ntracr if (trcflg(ktr).eq.0) then !100% in the mixed layer if (prij(k).le.dpmixl(i,j,1)) then tracer(i,j,k,1,ktr)=10.0 tracer(i,j,k,2,ktr)=10.0 else tracer(i,j,k,1,ktr)=0.0 tracer(i,j,k,2,ktr)=0.0 endif elseif (trcflg(ktr).eq.1) then !20 below euphotic zone if (swfrac(k).gt.0.01) then tracer(i,j,k,1,ktr)=0.0 tracer(i,j,k,2,ktr)=0.0 else tracer(i,j,k,1,ktr)=20.0 ! mg/m^3 tracer(i,j,k,2,ktr)=20.0 ! mg/m^3 endif elseif (trcflg(ktr).eq.2) then !temperature tracer(i,j,k,1,ktr)=temp(i,j,k,1) tracer(i,j,k,2,ktr)=temp(i,j,k,1) elseif (trcflg(ktr).eq.3) then !fully passive tracer(i,j,k,1,ktr)=0.0 !should never get here tracer(i,j,k,2,ktr)=0.0 !should never get here !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! inicializacao do NPDZ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! elseif (trcflg(ktr).eq.904 .or. & trcflg(ktr).eq.903 ) then !NPZD or NPZ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!! BIO_P !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! zk = 0.5*(prij(k+1)+prij(k))*qonem if (zk.le.300.0) then ! 0.1 at 300m, 1.0 at 100m, 2.025 at 0m ! bio_p = 0.1 + (300.0-zk)**2 * (0.9/200.0**2) bio_p = 0.01 + (300.0-zk)**2 * (0.025/200.0**2) !erro1 elseif (zk.le.900.0) then ! 0.1 at 300m, 0.0 at 900m ! bio_p = (900.0-zk) * 0.1/600.0 !erro2 bio_p = (900.0-zk) * 0.01/600.0 else bio_p = 0.0 endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!! BIO_N !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! if (temp(i,j,k,1).lt. 6.0) then ! bio_n = 37.0 ! elseif (temp(i,j,k,1).gt.27.0) then ! bio_n = 0.0 ! else ! bio_n = (27.0-temp(i,j,k,1)) * 37.0/21.0 ! bio_n = 39.3116-1.335*temp(i,j,k,1) ! endif !!!! variação linear segundo climatologia WOA13 para domínio do 1/24 !!!! (Daniel e Afonso - 7 de julho de 2014) - Luana 10/12/2021 ! zk = 0.5*(p(i,j,k+1)+p(i,j,k))*qonem !erro3 if (zk.le.100.0) then bio_n = (0.3*(100.-zk) + 0.9*(zk-0.)) & /(100.-0.) elseif (zk.le.800.0) then bio_n = (0.9*(800.-zk) + 30.*(zk-100.)) & /(800.-100.) elseif (zk.le.1200.0) then bio_n = (30.*(1200.-zk) + 32.*(zk-800.)) & /(1200.-800.) elseif (zk.le.2200.0) then bio_n = (32.*(2200.-zk) + 21.*(zk-1200.)) & /(2200.-1200.) else bio_n = 21. endif !!! modelando o detrito com valores do NTP do HABITATS !!! 0.4 na superficie e 70% a menos abaixo dos 100 m ! if (zk.le.100.0) then ! 0.4-(0.7*0.4) at 100m, 0.4 at 0m ! bio_d = (0.4*(100.-zk) + (0.4-0.28)*(zk-0.)) ! & /(100.-0.) ! else ! bio_d = 0.4-0.28 ! bio_d = 0.1 ! endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! tracer(i,j,k,1,ktr )=bio_n !N tracer(i,j,k,2,ktr )=bio_n tracer(i,j,k,1,ktr+1)=bio_p !P tracer(i,j,k,2,ktr+1)=bio_p tracer(i,j,k,1,ktr+2)=bio_p !Z=P tracer(i,j,k,2,ktr+2)=bio_p if (trcflg(ktr).eq.904) then ! tracer(i,j,k,1,ktr+3)=bio_p + 1.0 !D=P+1 ! tracer(i,j,k,2,ktr+3)=bio_p + 1.0 tracer(i,j,k,1,ktr+3)=0.1 ! segundo sugestões diversas (Lima e outros) tracer(i,j,k,2,ktr+3)=0.1 endif endif !trcflg !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! enddo !ktr enddo !k endif !ip enddo !i enddo !j !$OMP END PARALLEL DO endif !iniflg.eq.2:else ! if (lpipe .and. lpipe_initrc) then do ktr= 1,ntracr do k= 1,kk write (ptxt,'(a4,i2.2,a3,i3)') 'trc.',ktr,' k=',k call pipe_compare_sym1(tracer(1-nbdy,1-nbdy,k,1,ktr), & ip,ptxt) enddo !k enddo !ktr endif !lpipe.and.lpipe_initrc ! if (itest.gt.0 .and. jtest.gt.0) then prij(1)=0.0 do k=1,kk prij(k+1)=prij(k)+dp(itest,jtest,k,1) enddo !k write(cformat,'(a,i2,a,i2,a)') & '(i9,2i5,a,',ntracr, & 'a / (23x,i3,2f8.2,', ntracr,'f8.4))' write (lp,cformat) & nstep,i0+itest,j0+jtest, & ' istate: thkns dpth', & (' tracer',ktr=1,ntracr), & (k, & dp(itest,jtest,k,1)*qonem, & (prij(k+1)+prij(k))*0.5*qonem, & (tracer(itest,jtest,k,1,ktr),ktr=1,ntracr), & k=1,kk) write(lp,'(23x,a,8x,f8.2)') 'bot',depths(itest,jtest) endif !test tile call xcsync(flush_lp) ! return end subroutine trcupd(m,n) use mod_xc ! HYCOM communication interface use mod_cb_arrays ! HYCOM saved arrays implicit none ! integer m,n ! ! --- ----------------------------------------------------------- ! --- tracer-specific operations (side-wall relaxation in thermf) ! --- ----------------------------------------------------------- ! integer i,j,k,ktr real chl,pij(kdm+1),swfrac(kdm+1),q ! do ktr= 1,ntracr if (trcflg(ktr).eq.0) then if (trcrlx) then ! --- tracer always trwall, when non-zero, at surface !$OMP PARALLEL DO PRIVATE(j,k,i,q) & !$OMP SCHEDULE(STATIC,jblk) do j=1,jj do i=1,ii if (SEA_P) then ! --- constant in time trwall remains exactly constant q = trwall(i,j,1,lc0,ktr) & +(trwall(i,j,1,lc1,ktr)-trwall(i,j,1,lc0,ktr))*wc1 & +(trwall(i,j,1,lc2,ktr)-trwall(i,j,1,lc0,ktr))*wc2 & +(trwall(i,j,1,lc3,ktr)-trwall(i,j,1,lc0,ktr))*wc3 if (q.gt.0.0) then tracer(i,j,1,n,ktr) = q endif endif !ip enddo !i enddo !j !$OMP END PARALLEL DO elseif (.not. trcrlx) then ! --- tracer always 10.0 at surface !$OMP PARALLEL DO PRIVATE(j,k,i) & !$OMP SCHEDULE(STATIC,jblk) do j=1,jj do i=1,ii if (SEA_P) then tracer(i,j,1,n,ktr) = 10.0 endif !ip enddo !i enddo !j !$OMP END PARALLEL DO endif !trcrlx:else elseif (trcflg(ktr).eq.1) then ! --- psudo-silicate, half-life of 30 days in euphotic zone q = 1.0-delt1/(30.0*86400.0) !$OMP PARALLEL DO PRIVATE(j,k,i,chl,pij,swfrac) & !$OMP SCHEDULE(STATIC,jblk) do j=1,jj do i=1,ii if (SEA_P) then pij(1)=0.0 do k=1,kk pij(k+1) = pij(k)+dp(i,j,k,n) enddo !k if (jerlv0.le.0) then chl = akpar(i,j,lk0)*wk0+akpar(i,j,lk1)*wk1 & +akpar(i,j,lk2)*wk2+akpar(i,j,lk3)*wk3 endif call swfrac_ij(chl,pij,kk+1,qonem*oneta(i,j,n), & jerlov(i,j),swfrac) do k=1,kk if (0.5*(swfrac(k)+swfrac(k+1)).gt.0.01) then tracer(i,j,k,n,ktr) = q*tracer(i,j,k,n,ktr) else exit !too deep endif enddo !k endif !ip enddo !i enddo !j !$OMP END PARALLEL DO elseif (trcflg(ktr).eq.2) then ! --- temperature-like (do nothing, heat flux forcing in mixed layer) elseif (trcflg(ktr).eq.3) then ! --- fully passive (do nothing) elseif (trcflg(ktr).eq.903) then ! --- NPZ call trcupd_903(m,n, ktr) elseif (trcflg(ktr).eq.904) then ! --- NPZD call trcupd_904(m,n, ktr) ! elseif (trcflg(ktr).eq.909) then ! --- Chai 9-component. ! call trcupd_909(m,n, ktr) endif enddo !ktr return end subroutine trcupd subroutine trcupd_903(m,n, ibio) use mod_xc ! HYCOM communication interface use mod_cb_arrays ! HYCOM saved arrays implicit none ! integer m,n,ibio ! ! --- ------------------------------------------------- ! --- tracer-specific operations for Franks NPZ biology ! --- ------------------------------------------------- ! real, save, dimension(mxtrcr) :: & bup, & ! maximum growth rate of phytoplankton (1/d). bgz, & ! maximum grazing rate of zooplankton (1/d). bdp, & ! senescence (death) rate of phytoplankton (1/d). bdz, & ! death rate of zooplankton (1/d). buk, & ! = half-saturation coefficient for phytoplankton (mg/m^3) asim, & ! assimilation efficiency of zooplankton. glam ! Ivlev parameter for grazing efficiency of zooplankton. ! integer i,j,k real bm_n,bm_p,bm_z,bn_n,bn_p,bn_z,bu_n,bu_p,bu_z, & uptake,grazin,pdeath,zdeath, & chl,par,pij(kdm+1),swfrac(kdm+1) ! if (ibio.lt.0) then !initialize only ! ! --- read from tracer_NN.input: ! --- 'biotyp' = type (90X=std.bio,X=3,4,9) must be 903 ! --- 'bup ' = maximum growth rate of phytoplankton (1/d). ! --- 'bgz ' = maximum grazing rate of zooplankton (1/d). ! --- 'bdp ' = senescence (death) rate of phytoplankton (1/d). ! --- 'bdz ' = death rate of zooplankton (1/d). ! --- 'buk ' = half-saturation coefficient for phytoplankton (mg/m^3) ! --- 'asim ' = assimilation efficiency of zooplankton. ! --- 'glam ' = Ivlev parameter for grazing efficiency of zooplankton. ! i = -ibio if (mnproc.eq.1) then write(lp,'(/ a,i3,a,i3,a)') & 'Franks NPZ parameters for tracers',i,' to',i+2,':' endif !1st tile ! call blkini(k, 'biotyp') if (k.ne.903) then if (mnproc.eq.1) then write(lp,'(/ a /)') & 'error - biotyp must be 903' call flush(lp) endif !1st tile call xcstop('(trcini)') stop '(trcini)' endif !biotyp.ne.903 ! call blkinr(bup( i), 'bup ','(a6," =",f10.4," 1/d")') call blkinr(bgz( i), 'bgz ','(a6," =",f10.4," 1/d")') call blkinr(bdp( i), 'bdp ','(a6," =",f10.4," 1/d")') call blkinr(bdz( i), 'bdz ','(a6," =",f10.4," 1/d")') call blkinr(buk( i), 'buk ','(a6," =",f10.4," mg/m^3")') call blkinr(asim( i), 'asim ','(a6," =",f10.4," ")') call blkinr(glam( i), 'glam ','(a6," =",f10.4," ")') ! if (mnproc.eq.1) then write(lp,*) endif !1st tile return endif !ibio.lt.0 ! ! --- leapfrog time step. ! !$OMP PARALLEL DO PRIVATE(j,i,k,chl,pij,par,swfrac, & !$OMP bm_n,bm_p,bm_z,bn_n,bn_p,bn_z, & !$OMP bu_n,bu_p,bu_z, & !$OMP uptake,grazin,pdeath,zdeath) & !$OMP SCHEDULE(STATIC,jblk) do j=1,jj do i=1,ii if (SEA_P) then pij(1)=0.0 do k=1,kk pij(k+1) = pij(k)+dp(i,j,k,n) enddo !k if (jerlv0.le.0) then chl = akpar(i,j,lk0)*wk0+akpar(i,j,lk1)*wk1 & +akpar(i,j,lk2)*wk2+akpar(i,j,lk3)*wk3 endif call swfrac_ij(chl,pij,kk+1,qonem*oneta(i,j,n), & jerlov(i,j),swfrac) do k=1,kk par = 0.5*(swfrac(k)+swfrac(k+1)) ! bm_n = tracer(i,j,k,m,ibio) bm_p = tracer(i,j,k,m,ibio+1) bm_z = tracer(i,j,k,m,ibio+2) bn_n = tracer(i,j,k,n,ibio) bn_p = tracer(i,j,k,n,ibio+1) bn_z = tracer(i,j,k,n,ibio+2) ! uptake = bup(ibio)*bm_p*bm_n*par/(buk(ibio)+bm_n) grazin = bgz(ibio)*bm_z*(1.0-exp(-glam(ibio)*bm_p)) pdeath = bdp(ibio)*bm_p zdeath = bdz(ibio)*bm_z ! limit negative terms to 10% of total per single time step grazin = min(grazin,bn_p*0.2*86400.0/delt1) uptake = min(uptake,bn_n*0.2*86400.0/delt1) ! bu_p = -grazin +uptake-pdeath bu_z = asim(ibio) *grazin-zdeath bu_n = (1.0-asim(ibio))*grazin+zdeath-uptake+pdeath ! tracer(i,j,k,n,ibio) = bn_n + delt1/86400.0 * bu_n tracer(i,j,k,n,ibio+1) = bn_p + delt1/86400.0 * bu_p tracer(i,j,k,n,ibio+2) = bn_z + delt1/86400.0 * bu_z ! ! --- fields must be non-negative ! --- note: only round-off should make a field negative ! if (tracer(i,j,k,n,ibio+1).lt.0.0) then !PtoN tracer(i,j,k,n,ibio) = tracer(i,j,k,n,ibio) - & tracer(i,j,k,n,ibio+1) tracer(i,j,k,n,ibio+1) = 0.0 endif if (tracer(i,j,k,n,ibio+2).lt.0.0) then !ZtoN tracer(i,j,k,n,ibio) = tracer(i,j,k,n,ibio) - & tracer(i,j,k,n,ibio+2) tracer(i,j,k,n,ibio+2) = 0.0 endif if (tracer(i,j,k,n,ibio) .lt.0.0) then !NtoPZ (do last) tracer(i,j,k,n,ibio+1) = tracer(i,j,k,n,ibio+1) - & tracer(i,j,k,n,ibio)*0.5 tracer(i,j,k,n,ibio+2) = tracer(i,j,k,n,ibio+2) - & tracer(i,j,k,n,ibio)*0.5 tracer(i,j,k,n,ibio) = 0.0 endif enddo !k endif !ip enddo !i enddo !j !$OMP END PARALLEL DO return end subroutine trcupd_903 subroutine trcupd_904(m,n, ibio) use mod_xc ! HYCOM communication interface use mod_cb_arrays ! HYCOM saved arrays implicit none ! integer m,n,ibio ! ! --- ------------------------------------------------------- ! --- tracer-specific operations for Lima/Idrisi NPZD biology ! --- ------------------------------------------------------- ! real, save, dimension(mxtrcr) :: & pp, & ! zoopl: preference term for phytoplankton pz, & ! zoopl: preference term for zooplankton pd, & ! zoopl: preference term for detritus aa, & ! zoopl: assimilation efficiency am, & ! zoopl: metabolic efficiency fkz, & ! zoopl: half-saturation coefficient (mg/m^3) gmax, & ! zoopl: maximum growth rate (1/day) zmor ! zoopl: mortality (1/day) ! real, save, dimension(mxtrcr) :: & ! ik, & ! phyto: light absorption efficiency scalar (einst/m^2/h) fkp, & ! phyto: half-saturation coefficient (mg/m^3) pmax, & ! phyto: maximum growth rate (1/day) psen ! phyto: senescence (1/day) ! real, save, dimension(mxtrcr) :: & remn ! detri: remineralization (1/day) ! integer, save, dimension(mxtrcr) :: & spcflg ! tmpfn: species type (0=none,1=cold-water,2=warm-water) ! real, parameter :: & ! temperature function for cold-water species ! thornton and lessem (1978) theta1 = 16.0, & ! dependence on lower optimum temperature curve theta2 = 9.0, & ! dependence on higher optimum temperature curve theta3 = 11.0, & ! maximum temperature (upper tolerance level) q10l = 2.0, & ! the metabolic q10 for temperature response xk1 = 0.5, & ! scalar constant xk2 = 0.98, & ! scalar constant xk3 = 0.01, & ! scalar constant xk4 = 0.01 ! scalar constant ! real, parameter :: & ! temperature function for warm-water species tmax = 27.0, & ! Tfunc: maximum tolerated temperature topt = 25.0, & ! Tfunc: optimum temperature q10w = 2.0 ! Tfunc: the metabolic q10 for temperature response ! integer i,j,k real bm_n,bm_p,bm_z,bm_d,bn_n,bn_p,bn_z,bn_d, & bu_n,bu_p,bu_z,bu_d, & gamma1,gamma2,xnum,xkatheta,ynum,xkbtheta, & tijk,tfn,vw,xw,yw,zw, & pgrw,zgrw,pref,prf2,qprf,ztgx,dofz,pofz,zofz, & chl,par,pij(kdm+1),swfrac(kdm+1) ! if (ibio.lt.0) then !initialize only ! ! --- read from tracer.input: ! --- 'biotyp' = type (90X=std.bio,X=3,4,9) must be 904 ! ! --- 'pp ' = zoopl: preference term for phytoplankton ! --- 'pz ' = zoopl: preference term for zooplankton ! --- 'pd ' = zoopl: preference term for detritus ! --- 'aa ' = zoopl: assimilation efficiency ! --- 'am ' = zoopl: metabolic efficiency ! --- 'fkz ' = zoopl: half-saturation coefficient (mg/m^3) ! --- 'gmax ' = zoopl: maximum growth rate (1/day) ! --- 'zmor ' = zoopl: mortality (1/day) ! ! --- 'ik ' = phyto: light absorption efficiency scalar (einst/m^2/h) ! --- 'fkp ' = phyto: half-saturation coefficient (mg/m^3) ! --- 'pmax ' = phyto: maximum growth rate (1/day) ! --- 'psen ' = phyto: senescence (1/day) ! ! --- 'remn ' = detri: remineralization (1/day) ! ! --- 'spcflg' = tmpfn: species type (0=none,1=cold-water,2=warm-water) ! i = -ibio if (mnproc.eq.1) then write(lp,'(/ a,i3,a,i3,a)') & 'Lima/Idrisi NPZD parameters for tracers',i,' to',i+3,':' endif !1st tile ! call blkini(k, 'biotyp') if (k.ne.904) then if (mnproc.eq.1) then write(lp,'(/ a /)') & 'error - biotyp must be 904' call flush(lp) endif !1st tile call xcstop('(trcini)') stop '(trcini)' endif !biotyp.ne.904 ! call blkinr(pp( i), 'pp ','(a6," =",f10.4," ")') call blkinr(pz( i), 'pz ','(a6," =",f10.4," ")') call blkinr(pd( i), 'pd ','(a6," =",f10.4," ")') call blkinr(aa( i), 'aa ','(a6," =",f10.4," ")') call blkinr(am( i), 'am ','(a6," =",f10.4," ")') call blkinr(fkz( i), 'fkz ','(a6," =",f10.4," mg/m^3")') call blkinr(gmax( i), 'gmax ','(a6," =",f10.4," 1/day")') call blkinr(zmor( i), 'zmor ','(a6," =",f10.4," 1/day")') ! call blkinr(fkp( i), 'fkp ','(a6," =",f10.4," mg/m^3")') call blkinr(pmax( i), 'pmax ','(a6," =",f10.4," 1/day")') call blkinr(psen( i), 'psen ','(a6," =",f10.4," 1/day")') ! call blkinr(remn( i), 'remn ','(a6," =",f10.4," 1/day")') ! call blkini(spcflg(i),'spcflg') ! if (mnproc.eq.1) then write(lp,*) endif !1st tile return endif !ibio.lt.0 ! ! --- leapfrog time step. ! !$OMP PARALLEL DO PRIVATE(j,i,k,chl,pij,par,swfrac, & !$OMP bm_n,bm_p,bm_z,bm_d,bn_n,bn_p,bn_z,bn_d, & !$OMP bu_n,bu_p,bu_z,bu_d, & !$OMP gamma1,gamma2,xnum,xkatheta,ynum,xkbtheta, & !$OMP tijk,tfn,vw,xw,yw,zw, & !$OMP pgrw,zgrw,pref,prf2,qprf,ztgx,dofz,pofz,zofz) & !$OMP SCHEDULE(STATIC,jblk) do j=1,jj do i=1,ii if (SEA_P) then pij(1)=0.0 do k=1,kk pij(k+1) = pij(k)+dp(i,j,k,n) enddo !k if (jerlv0.le.0) then chl = akpar(i,j,lk0)*wk0+akpar(i,j,lk1)*wk1 & +akpar(i,j,lk2)*wk2+akpar(i,j,lk3)*wk3 endif call swfrac_ij(chl,pij,kk+1,qonem*oneta(i,j,n), & jerlov(i,j),swfrac) do k=1,kk par = 0.5*(swfrac(k)+swfrac(k+1)) ! bm_n = tracer(i,j,k,m,ibio) bm_p = tracer(i,j,k,m,ibio+1) bm_z = tracer(i,j,k,m,ibio+2) bm_d = tracer(i,j,k,m,ibio+3) bn_n = tracer(i,j,k,n,ibio) bn_p = tracer(i,j,k,n,ibio+1) bn_z = tracer(i,j,k,n,ibio+2) bn_d = tracer(i,j,k,n,ibio+3) ! if (spcflg(ibio).eq.1) then ! --- cold-water species temperature dependance tijk = temp(i,j,k,n) gamma1 = 1.0/(theta2-q10l) * & log((xk2*(1.0-xk1))/(xk1*(1.0-xk2))) gamma2 = 1.0/(theta1-theta3) * & log((xk2*(1.0-xk3))/(xk4*(1.0-xk2))) xnum = exp(gamma1*(tijk-q10l)) xkatheta = (xk1*xnum)/(1.0+xk1*(xnum-1.0)) ynum = exp(gamma2*(theta1-tijk)) xkbtheta = (xk4*ynum)/(1.0+xk3*(ynum-1.0)) tfn = xkatheta*xkbtheta elseif (spcflg(ibio).eq.2) then ! --- warm-water species temperature dependance tijk = temp(i,j,k,n) if (tijk.le.tmax) then vw = (tmax-tijk)/(tmax-topt) yw = log(q10w)*(tmax-topt+2.0) zw = log(q10w)*(tmax-topt) xw = (zw**2 * (1.0+sqrt(1.0+40.0/yw))**2)/400.0 tfn = vw**xw * exp(xw*(1.0-vw)) else tfn=0.0 endif else ! --- no temperature dependance tfn=1.0 endif !spcflg ! pref = pp(ibio)*bm_p + & pd(ibio)*bm_d + & pz(ibio)*bm_z prf2 = pp(ibio)*bm_p**2 + & pd(ibio)*bm_d**2 + & pz(ibio)*bm_z**2 qprf = 1.0/(fkz(ibio)*pref + prf2 + epsil) !epsil prevents 1/0 ztgx = bm_z*tfn*gmax(ibio) ! pgrw = bm_p*tfn*pmax(ibio)*bm_n*par/(fkp(ibio)+bm_n) zgrw = ztgx*(prf2 *qprf)*aa(ibio)*am(ibio) pofz = ztgx*(pp(ibio)*bm_p**2*qprf) zofz = ztgx*(pz(ibio)*bm_z**2*qprf) dofz = ztgx*(pd(ibio)*bm_d**2*qprf) ! ! limit negative terms to 10% of total per single time step pgrw = min(pgrw,bn_n*0.2*86400.0/delt1) zgrw = min(zgrw,bn_n*0.2*86400.0/delt1) pofz = min(pofz,bn_p*0.2*86400.0/delt1) zofz = min(zofz,bn_z*0.2*86400.0/delt1) dofz = min(dofz,bn_d*0.2*86400.0/delt1) ! bu_p = pgrw & - pofz & - bm_p*psen(ibio) bu_z = zgrw & - zofz & - bm_z**2*zmor(ibio) ! decaimento quadratico daniel 23/10/14 !erro4 ! - bm_z*zmor(ibio) bu_d = bm_p*psen(ibio) & ! + bm_z*zmor(ibio) & ! + (pofz+zofz+dofz)*(1.0-aa(ibio)) & + (bm_z**2)*zmor(ibio) & ! decaimento quadratico !erro5 + (pofz+zofz+dofz)*(1.0-aa(ibio)) & - dofz & - bm_d*remn(ibio) bu_n = bm_d*remn(ibio) & + (pofz+zofz+dofz)* aa(ibio) & - zgrw & - pgrw ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!! afonso - Luana ! a varaiavel ibi passa no common para salver tracer biologico no mooring ibi = ibio !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! tracer(i,j,k,n,ibio) = bn_n + delt1/86400.0 * bu_n tracer(i,j,k,n,ibio+1) = bn_p + delt1/86400.0 * bu_p tracer(i,j,k,n,ibio+2) = bn_z + delt1/86400.0 * bu_z tracer(i,j,k,n,ibio+3) = bn_d + delt1/86400.0 * bu_d ! ! --- fields must be non-negative ! --- note: only round-off should make a field negative ! if (tracer(i,j,k,n,ibio+1).lt.0.0) then !PtoN tracer(i,j,k,n,ibio) = tracer(i,j,k,n,ibio) - & tracer(i,j,k,n,ibio+1) tracer(i,j,k,n,ibio+1) = 0.0 endif if (tracer(i,j,k,n,ibio+2).lt.0.0) then !ZtoN tracer(i,j,k,n,ibio) = tracer(i,j,k,n,ibio) - & tracer(i,j,k,n,ibio+2) tracer(i,j,k,n,ibio+2) = 0.0 endif if (tracer(i,j,k,n,ibio+3).lt.0.0) then !DtoN tracer(i,j,k,n,ibio) = tracer(i,j,k,n,ibio) - & tracer(i,j,k,n,ibio+3) tracer(i,j,k,n,ibio+3) = 0.0 endif if (tracer(i,j,k,n,ibio) .lt.0.0) then !NtoD (do last) tracer(i,j,k,n,ibio+3) = tracer(i,j,k,n,ibio+3) - & tracer(i,j,k,n,ibio) tracer(i,j,k,n,ibio) = 0.0 endif enddo endif !ip enddo !i enddo !j !$OMP END PARALLEL DO return end subroutine trcupd_904 subroutine pcmtrc(si,pi,ki,ks, so,po,ko) implicit none ! integer ki,ks,ko real si(ki,ks),pi(ki+1), & so(ko,ks),po(ko+1) ! !********** !* ! 1) remap from one set of vertical cells to another. ! method: piecewise constant across each input cell ! the output is the average of the interpolation ! profile across each output cell. ! ! 2) input arguments: ! si - scalar fields in pi-layer space ! pi - layer interface depths (non-negative m) ! pi( 1) is the surface ! pi(ki+1) is the bathymetry ! ki - 1st dimension of si (number of input layers) ! ks - 2nd dimension of si,so (number of fields) ! po - target interface depths (non-negative m) ! po(k+1) >= po(k) ! ko - 1st dimension of so (number of output layers) ! ! 3) output arguments: ! so - scalar fields in po-layer space ! ! 4) except at data voids, must have: ! pi( 1) == zero (surface) ! pi( l+1) >= pi(l) ! pi(ki+1) == bathymetry ! 0 <= po(k) <= po(k+1) ! output layers completely below the bathymetry inherit values ! from the layer above. ! ! 5) Alan J. Wallcraft, Naval Research Laboratory, Sep. 2002 (Aug. 2005). !* !********** ! real thin parameter (thin=1.e-6) ! minimum layer thickness (no division by 0.0) ! integer i,k,l,lf real q,zb,zt,sok(ks) ! lf=1 zb=po(1) do k= 1,ko zt = zb zb = po(k+1) ! WRITE(6,*) 'k,zt,zb = ',k,zt,zb if (zb-zt.lt.thin .or. zt.ge.pi(ki+1)) then ! ! --- thin or bottomed layer, values taken from layer above ! do i= 1,ks so(k,i) = so(k-1,i) enddo !i else ! ! form layer averages. ! if (pi(lf).gt.zt) then WRITE(6,*) 'bad lf = ',lf stop endif do i= 1,ks sok(i) = 0.0 enddo !i do l= lf,ki if (pi(l).gt.zb) then ! WRITE(6,*) 'l,lf= ',l,lf,l-1 lf = l-1 exit elseif (pi(l).ge.zt .and. pi(l+1).le.zb) then ! ! the input layer is completely inside the output layer ! q = max(pi(l+1)-pi(l),thin)/(zb-zt) do i= 1,ks sok(i) = sok(i) + q*si(l,i) enddo !i ! WRITE(6,*) 'L,q = ',l,q else ! ! the input layer is partially inside the output layer ! q = max(min(pi(l+1),zb)-max(pi(l),zt),thin)/(zb-zt) do i= 1,ks sok(i) = sok(i) + q*si(l,i) enddo !i ! WRITE(6,*) 'l,q = ',l,q endif enddo !l do i= 1,ks so(k,i) = sok(i) enddo !i endif enddo !k return end subroutine pcmtrc subroutine plmtrc(si,pi,ki,ks, so,po,ko) implicit none ! integer ki,ks,ko real si(ki,ks),pi(ki+1), & so(ko,ks),po(ko+1),flag ! !********** !* ! 1) remap from one set of vertical cells to another. ! method: piecewise linear across each input cell ! the output is the average of the interpolation ! profile across each output cell. ! ! 2) input arguments: ! si - scalar fields in pi-layer space ! pi - layer interface depths (non-negative m) ! pi( 1) is the surface ! pi(ki+1) is the bathymetry ! ki - 1st dimension of si (number of input layers) ! ks - 2nd dimension of si,so (number of fields) ! po - target interface depths (non-negative m) ! po(k+1) >= po(k) ! ko - 1st dimension of so (number of output layers) ! flag - data void (land) marker ! ! 3) output arguments: ! so - scalar fields in po-layer space ! ! 4) except at data voids, must have: ! pi( 1) == zero (surface) ! pi( l+1) >= pi(l) ! pi(ki+1) == bathymetry ! 0 <= po(k) <= po(k+1) ! output layers completely below the bathymetry inherit values ! from the layer above. ! ! 5) Tim Campbell, Mississippi State University, October 2002. ! Alan J. Wallcraft, Naval Research Laboratory, Aug. 2005. !* !********** ! real,parameter :: thin=1.e-6 !minimum layer thickness ! integer i,k,l,lf real q,qc,zb,zc,zt,sok(ks) real sis(ki,ks),pit(ki+1) real si_min(ks),si_max(ks) ! ! --- inforce minval(si(:,i)) <= minval(so(:,i)) and ! --- maxval(si(:,i)) >= maxval(so(:,i)) for i=1:ks ! --- in particular this inforces non-negativity, e.g. of tracers ! --- only required due to finite precision ! do i= 1,ks si_min(i) = minval(si(:,i)) si_max(i) = maxval(si(:,i)) enddo !i ! ! --- compute PLM slopes for input layers do k=1,ki pit(k)=max(pi(k+1)-pi(k),thin) enddo call plmtrcx(pit,si,sis,ki,ks) ! --- compute output layer averages lf=1 zb=po(1) do k= 1,ko zt = zb zb = po(k+1) ! WRITE(6,*) 'k,zt,zb = ',k,zt,zb if (zb-zt.lt.thin .or. zt.ge.pi(ki+1)) then ! ! --- thin or bottomed layer, values taken from layer above ! do i= 1,ks so(k,i) = so(k-1,i) enddo !i else ! ! form layer averages. ! if (pi(lf).gt.zt) then WRITE(6,*) 'bad lf = ',lf stop endif do i= 1,ks sok(i) = 0.0 enddo !i do l= lf,ki if (pi(l).gt.zb) then ! WRITE(6,*) 'l,lf= ',l,lf,l-1 lf = l-1 exit elseif (pi(l).ge.zt .and. pi(l+1).le.zb) then ! ! the input layer is completely inside the output layer ! q = max(pi(l+1)-pi(l),thin)/(zb-zt) do i= 1,ks sok(i) = sok(i) + q*si(l,i) enddo !i ! WRITE(6,*) 'L,q = ',l,q else ! ! the input layer is partially inside the output layer ! average of linear profile is its center value ! q = max( min(pi(l+1),zb)-max(pi(l),zt), thin )/(zb-zt) zc = 0.5*(min(pi(l+1),zb)+max(pi(l),zt)) qc = (zc-pi(l))/pit(l) - 0.5 do i= 1,ks sok(i) = sok(i) + q*(si(l,i) + qc*sis(l,i)) enddo !i ! WRITE(6,*) 'l,q,qc = ',l,q,qc endif enddo !l do i= 1,ks so(k,i) = sok(i) so(k,i) = max( si_min(i), so(k,i) ) so(k,i) = min( si_max(i), so(k,i) ) enddo !i endif enddo !k return end subroutine plmtrc subroutine plmtrcx(pt, s,ss,ki,ks) implicit none ! integer ki,ks real pt(ki+1),s(ki,ks),ss(ki,ks) ! !********** !* ! 1) generate a monotonic PLM interpolation of a layered field ! ! 2) input arguments: ! pt - layer interface thicknesses (non-zero) ! s - scalar fields in layer space ! ki - 1st dimension of s (number of layers) ! ks - 2nd dimension of s (number of fields) ! ! 3) output arguments: ! ss - scalar field slopes for PLM interpolation ! ! 4) except at data voids, must have: ! pi( 1) == zero (surface) ! pi( l+1) >= pi(:,:,l) ! pi(ki+1) == bathymetry ! ! 5) Tim Campbell, Mississippi State University, September 2002. !* !********** ! integer l real ql(ki),qc(ki),qr(ki) ! !compute grid spacing ratios for slope computations ql(1)=0.0 qc(1)=0.0 qr(1)=0.0 do l=2,ki-1 ql(l)=2.0*pt(l)/(pt(l-1)+pt(l)) qc(l)=2.0*pt(l)/(pt(l-1)+2.0*pt(l)+pt(l+1)) qr(l)=2.0*pt(l)/(pt(l)+pt(l+1)) enddo ql(ki)=0.0 qc(ki)=0.0 qr(ki)=0.0 !compute normalized layer slopes do l=1,ks call plmtrcs(ql,qc,qr,s(1,l),ss(1,l),ki) enddo return end subroutine plmtrcx subroutine plmtrcs(rl,rc,rr,a,s,n) implicit none ! integer,intent(in) :: n real, intent(in) :: rl(n),rc(n),rr(n),a(n) real, intent(out) :: s(n) ! !********** !* ! 1) generate slopes for monotonic piecewise linear distribution ! ! 2) input arguments: ! rl - left grid spacing ratio ! rc - center grid spacing ratio ! rr - right grid spacing ratio ! a - scalar field zone averages ! n - number of zones ! ! 3) output arguments: ! s - zone slopes ! ! 4) Tim Campbell, Mississippi State University, September 2002. !* !********** ! integer,parameter :: ic=2, im=1, imax=100 real,parameter :: fracmin=1e-6, dfac=0.5 ! integer i,j real sl,sc,sr real dnp,dnn,dl,dr,ds,frac ! ! Compute zone slopes ! Campbell Eq(15) -- nonuniform grid ! s(1)=0.0 do j=2,n-1 sl=rl(j)*(a(j)-a(j-1)) sr=rr(j)*(a(j+1)-a(j)) if (sl*sr.gt.0.) then s(j)=sign(min(abs(sl),abs(sr)),sl) else s(j)=0.0 endif enddo s(n)=0.0 ! ! Minimize discontinuities between zones ! Apply single pass discontinuity minimization: Campbell Eq(19) ! do j=2,n-1 if(s(j).ne.0.0) then dl=-0.5*(s(j)+s(j-1))+a(j)-a(j-1) dr=-0.5*(s(j+1)+s(j))+a(j+1)-a(j) ds=sign(min(abs(dl),abs(dr)),dl) s(j)=s(j)+2.0*ds endif enddo return end subroutine plmtrcs ! ! !> Revision history: !> !> Aug 2002 - new routine to put all tracer interactions in one place !> Dec. 2003 - inforce non-negative bio-tracers !> Aug. 2005 - interpolate trwall to actual layer structure !> Aug. 2012 - constant in time trwall remains exactly constant !> Oct. 2013 - added jerlv0=-1 and calls to swfrac_ij !> May 2014 - use land/sea masks (e.g. ip) to skip land !> Nov. 2018 - allow for oneta in swfrac except in initrc