diff --git a/hamocc/mo_aufr_bgc.F90 b/hamocc/mo_aufr_bgc.F90 index 513592c2..8707ea62 100644 --- a/hamocc/mo_aufr_bgc.F90 +++ b/hamocc/mo_aufr_bgc.F90 @@ -95,7 +95,7 @@ subroutine aufr_bgc(kpie,kpje,kpke,ntr,ntrbgc,itrbgc,trc,kplyear,kplmon,kplday,o use mo_vgrid, only: kbo use mo_sedmnt, only: sedhpl use mo_intfcblom, only: sedlay2,powtra2,burial2,atm2 - use mo_param_bgc, only: bifr13,bifr14,c14fac,re1312,re14to,prei13,prei14 + use mo_param_bgc, only: bifr13_ini,bifr14_ini,c14fac,re1312,re14to,prei13,prei14 use mo_netcdf_bgcrw, only: read_netcdf_var ! Arguments @@ -518,14 +518,14 @@ subroutine aufr_bgc(kpie,kpje,kpke,ntr,ntrbgc,itrbgc,trc,kplyear,kplmon,kplday,o ! Initialise the remaining 13C and 14C fields, using the restart isco212 field rco213=locetra(i,j,k,isco213)/(locetra(i,j,k,isco212)+safediv) rco214=locetra(i,j,k,isco214)/(locetra(i,j,k,isco212)+safediv) - locetra(i,j,k,idoc13)=locetra(i,j,k,idoc)*rco213*bifr13 - locetra(i,j,k,idoc14)=locetra(i,j,k,idoc)*rco214*bifr14 - locetra(i,j,k,iphy13)=locetra(i,j,k,iphy)*rco213*bifr13 - locetra(i,j,k,iphy14)=locetra(i,j,k,iphy)*rco214*bifr14 - locetra(i,j,k,izoo13)=locetra(i,j,k,izoo)*rco213*bifr13 - locetra(i,j,k,izoo14)=locetra(i,j,k,izoo)*rco214*bifr14 - locetra(i,j,k,idet13)=locetra(i,j,k,idet)*rco213*bifr13 - locetra(i,j,k,idet14)=locetra(i,j,k,idet)*rco214*bifr14 + locetra(i,j,k,idoc13)=locetra(i,j,k,idoc)*rco213*bifr13_ini + locetra(i,j,k,idoc14)=locetra(i,j,k,idoc)*rco214*bifr14_ini + locetra(i,j,k,iphy13)=locetra(i,j,k,iphy)*rco213*bifr13_ini + locetra(i,j,k,iphy14)=locetra(i,j,k,iphy)*rco214*bifr14_ini + locetra(i,j,k,izoo13)=locetra(i,j,k,izoo)*rco213*bifr13_ini + locetra(i,j,k,izoo14)=locetra(i,j,k,izoo)*rco214*bifr14_ini + locetra(i,j,k,idet13)=locetra(i,j,k,idet)*rco213*bifr13_ini + locetra(i,j,k,idet14)=locetra(i,j,k,idet)*rco214*bifr14_ini locetra(i,j,k,icalc13)=locetra(i,j,k,icalc)*rco213 locetra(i,j,k,icalc14)=locetra(i,j,k,icalc)*rco214 endif @@ -542,8 +542,8 @@ subroutine aufr_bgc(kpie,kpje,kpke,ntr,ntrbgc,itrbgc,trc,kplyear,kplmon,kplday,o rco214=locetra(i,j,kbo(i,j),isco214)/(locetra(i,j,kbo(i,j),isco212)+safediv) powtra2(i,j,k,ipowc13)=powtra2(i,j,k,ipowaic)*rco213 powtra2(i,j,k,ipowc14)=powtra2(i,j,k,ipowaic)*rco214 - sedlay2(i,j,k,issso13)=sedlay2(i,j,k,issso12)*rco213*bifr13 - sedlay2(i,j,k,issso14)=sedlay2(i,j,k,issso12)*rco214*bifr14 + sedlay2(i,j,k,issso13)=sedlay2(i,j,k,issso12)*rco213*bifr13_ini + sedlay2(i,j,k,issso14)=sedlay2(i,j,k,issso12)*rco214*bifr14_ini sedlay2(i,j,k,isssc13)=sedlay2(i,j,k,isssc12)*rco213 sedlay2(i,j,k,isssc14)=sedlay2(i,j,k,isssc12)*rco214 endif @@ -557,8 +557,8 @@ subroutine aufr_bgc(kpie,kpje,kpke,ntr,ntrbgc,itrbgc,trc,kplyear,kplmon,kplday,o if(omask(i,j) > 0.5) then rco213=locetra(i,j,kbo(i,j),isco213)/(locetra(i,j,kbo(i,j),isco212)+safediv) rco214=locetra(i,j,kbo(i,j),isco214)/(locetra(i,j,kbo(i,j),isco212)+safediv) - burial2(i,j,k,issso13)=burial2(i,j,k,issso12)*rco213*bifr13 - burial2(i,j,k,issso14)=burial2(i,j,k,issso12)*rco214*bifr14 + burial2(i,j,k,issso13)=burial2(i,j,k,issso12)*rco213*bifr13_ini + burial2(i,j,k,issso14)=burial2(i,j,k,issso12)*rco214*bifr14_ini burial2(i,j,k,isssc13)=burial2(i,j,k,isssc12)*rco213 burial2(i,j,k,isssc14)=burial2(i,j,k,isssc12)*rco214 endif diff --git a/hamocc/mo_biomod.F90 b/hamocc/mo_biomod.F90 index ab381cc3..5ea92b33 100644 --- a/hamocc/mo_biomod.F90 +++ b/hamocc/mo_biomod.F90 @@ -84,9 +84,6 @@ module mo_biomod real, dimension (:,:), allocatable, public :: int_chbr3_prod real, dimension (:,:), allocatable, public :: int_chbr3_uv - real, public :: growth_co2 - real, public :: bifr13_perm - CONTAINS subroutine alloc_mem_biomod(kpie,kpje,kpke) diff --git a/hamocc/mo_carchm.F90 b/hamocc/mo_carchm.F90 index 90151b0d..c7fdb6f7 100644 --- a/hamocc/mo_carchm.F90 +++ b/hamocc/mo_carchm.F90 @@ -170,13 +170,13 @@ subroutine carchm(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,prho,pglat,omask,psicomo !$OMP ,kwco2,kwdms,kwo2,atco2,ato2,atn2,fluxd,fluxu,oxflux,tc_sat & !$OMP ,niflux,n2oflux,dmsflux,omega,supsat,undsa,dissol & !$OMP ,sch_11,sch_12,sch_sf,kw_11,kw_12,kw_sf,a_11,a_12,a_sf,flx11 & - !$OMP ,flx12,flxsf,atm_cfc11,atm_cfc12,atm_sf6 & + !$OMP ,flx12,flxsf,atm_cfc11,atm_cfc12,atm_sf6,fact & !$OMP ,natcu,natcb,natcc,natpco2,natfluxd,natfluxu,natomega & !$OMP ,natsupsat,natundsa,natdissol & !$OMP ,atco213,atco214,rco213,rco214,pco213,pco214,frac_aqg & !$OMP ,frac_dicg,flux13d,flux13u,flux14d,flux14u,dissol13,dissol14 & !$OMP ,flx_bromo,sch_bromo,kw_bromo,a_bromo,atbrf,Kb1,lsub & - !$OMP ,j,i) + !$OMP ,k,j,i,rrho,scn2,scn2o,kwn2,kwn2o) do k=1,kpke do j=1,kpje do i=1,kpie @@ -572,7 +572,7 @@ subroutine carchm(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,prho,pglat,omask,psicomo ! C14 decay in the sediment (could be moved to sediment part) if (use_cisonew .and. .not. use_sedbypass) then do k=1,ks - !$OMP PARALLEL DO PRIVATE(i) + !$OMP PARALLEL DO PRIVATE(i,j) do j=1,kpje do i=1,kpie if(omask(i,j).gt.0.5) then diff --git a/hamocc/mo_dipowa.F90 b/hamocc/mo_dipowa.F90 index 449260dd..85e94665 100644 --- a/hamocc/mo_dipowa.F90 +++ b/hamocc/mo_dipowa.F90 @@ -66,8 +66,7 @@ subroutine dipowa(kpie,kpje,kpke,omask,lspin) real :: tredsy(kpie,0:kpke,3) ! redsy for 'reduced system'? real :: aprior ! start value of oceanic tracer in bottom layer - !$OMP PARALLEL DO & - !$OMP&PRIVATE(i,k,iv,l,tredsy,sedb1,aprior,iv_oc) + !$OMP PARALLEL DO PRIVATE(i,k,iv,l,tredsy,sedb1,aprior,iv_oc) j_loop: do j=1,kpje k = 0 @@ -182,6 +181,7 @@ subroutine dipowa(kpie,kpje,kpke,omask,lspin) endif ! .not. lspin enddo j_loop + !$OMP END PARALLEL DO end subroutine dipowa diff --git a/hamocc/mo_ini_fields.F90 b/hamocc/mo_ini_fields.F90 index f4c2898d..5b7e48ce 100644 --- a/hamocc/mo_ini_fields.F90 +++ b/hamocc/mo_ini_fields.F90 @@ -87,7 +87,7 @@ subroutine ini_fields_ocean(kpaufr,kpie,kpje,kpke,kbnd,pddpo,prho,omask,pglon,pg !*********************************************************************************************** use mo_carbch, only: co2star,co3,hi,ocetra - use mo_param_bgc, only: fesoly,cellmass,fractdim,bifr13,bifr14,c14fac,re1312,re14to + use mo_param_bgc, only: fesoly,cellmass,fractdim,bifr13_ini,bifr14_ini,c14fac,re1312,re14to use mo_biomod, only: abs_oce use mo_control_bgc, only: rmasks,use_FB_BGC_OCE,use_cisonew,use_AGG,use_CFC,use_natDIC, & use_BROMO, use_sedbypass @@ -213,14 +213,14 @@ subroutine ini_fields_ocean(kpaufr,kpie,kpje,kpke,kbnd,pddpo,prho,omask,pglon,pg if (use_cisonew) then rco213=ocetra(i,j,k,isco213)/(ocetra(i,j,k,isco212)+safediv) rco214=ocetra(i,j,k,isco214)/(ocetra(i,j,k,isco212)+safediv) - ocetra(i,j,k,iphy13) =ocetra(i,j,k,iphy)*rco213*bifr13 - ocetra(i,j,k,iphy14) =ocetra(i,j,k,iphy)*rco214*bifr14 - ocetra(i,j,k,izoo13) =ocetra(i,j,k,izoo)*rco213*bifr13 - ocetra(i,j,k,izoo14) =ocetra(i,j,k,izoo)*rco214*bifr14 - ocetra(i,j,k,idoc13) =ocetra(i,j,k,idoc)*rco213*bifr13 - ocetra(i,j,k,idoc14) =ocetra(i,j,k,idoc)*rco214*bifr14 - ocetra(i,j,k,idet13) =ocetra(i,j,k,idet)*rco213*bifr13 - ocetra(i,j,k,idet14) =ocetra(i,j,k,idet)*rco214*bifr14 + ocetra(i,j,k,iphy13) =ocetra(i,j,k,iphy)*rco213*bifr13_ini + ocetra(i,j,k,iphy14) =ocetra(i,j,k,iphy)*rco214*bifr14_ini + ocetra(i,j,k,izoo13) =ocetra(i,j,k,izoo)*rco213*bifr13_ini + ocetra(i,j,k,izoo14) =ocetra(i,j,k,izoo)*rco214*bifr14_ini + ocetra(i,j,k,idoc13) =ocetra(i,j,k,idoc)*rco213*bifr13_ini + ocetra(i,j,k,idoc14) =ocetra(i,j,k,idoc)*rco214*bifr14_ini + ocetra(i,j,k,idet13) =ocetra(i,j,k,idet)*rco213*bifr13_ini + ocetra(i,j,k,idet14) =ocetra(i,j,k,idet)*rco214*bifr14_ini ocetra(i,j,k,icalc13)=ocetra(i,j,k,icalc)*rco213 ocetra(i,j,k,icalc14)=ocetra(i,j,k,icalc)*rco214 endif @@ -268,10 +268,10 @@ subroutine ini_fields_ocean(kpaufr,kpie,kpje,kpke,kbnd,pddpo,prho,omask,pglon,pg if (use_cisonew) then rco213=ocetra(i,j,kbo(i,j),isco213)/(ocetra(i,j,kbo(i,j),isco212)+safediv) rco214=ocetra(i,j,kbo(i,j),isco214)/(ocetra(i,j,kbo(i,j),isco212)+safediv) - powtra(i,j,k,ipowc13)=powtra(i,j,k,ipowaic)*rco213*bifr13 - powtra(i,j,k,ipowc14)=powtra(i,j,k,ipowaic)*rco214*bifr14 - sedlay(i,j,k,issso13)=sedlay(i,j,k,issso12)*rco213*bifr13 - sedlay(i,j,k,issso14)=sedlay(i,j,k,issso12)*rco214*bifr14 + powtra(i,j,k,ipowc13)=powtra(i,j,k,ipowaic)*rco213*bifr13_ini + powtra(i,j,k,ipowc14)=powtra(i,j,k,ipowaic)*rco214*bifr14_ini + sedlay(i,j,k,issso13)=sedlay(i,j,k,issso12)*rco213*bifr13_ini + sedlay(i,j,k,issso14)=sedlay(i,j,k,issso12)*rco214*bifr14_ini sedlay(i,j,k,isssc13)=sedlay(i,j,k,isssc12)*rco213 sedlay(i,j,k,isssc14)=sedlay(i,j,k,isssc12)*rco214 endif diff --git a/hamocc/mo_intfcblom.F90 b/hamocc/mo_intfcblom.F90 index dece38d3..58bdbc35 100644 --- a/hamocc/mo_intfcblom.F90 +++ b/hamocc/mo_intfcblom.F90 @@ -224,7 +224,6 @@ subroutine blom2hamocc(m,n,mm,nn) ! --- calculate pressure at interfaces (necesarry since p has ! --- not been calculated at restart) - !$OMP PARALLEL DO PRIVATE(k,kn,l,i) do k=1,kk kn=k+nn do j=1,jj @@ -235,7 +234,6 @@ subroutine blom2hamocc(m,n,mm,nn) enddo enddo enddo - !$OMP END PARALLEL DO ! --- ------------------------------------------------------------------ ! --- 2D fields @@ -334,13 +332,22 @@ subroutine blom2hamocc(m,n,mm,nn) do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) sedlay(i,j,k,:) = sedlay2(i,j,kn,:) powtra(i,j,k,:) = powtra2(i,j,kn,:) - burial(i,j,:) = burial2(i,j,n,:) enddo enddo enddo enddo !$OMP END PARALLEL DO + !$OMP PARALLEL DO PRIVATE(l,i) + do j=1,jj + do l=1,isp(j) + do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) + burial(i,j,:) = burial2(i,j,n,:) + enddo + enddo + enddo + !$OMP END PARALLEL DO + endif ! --- ------------------------------------------------------------------ @@ -436,15 +443,25 @@ subroutine hamocc2blom(m,n,mm,nn) powtra2(i,j,km,:) = wts1*powtra2(i,j,km,:) & & + wts2*powtra2(i,j,kn,:) & & + wts2*powtra(i,j,k,:) - burial2(i,j,m,:) = wts1*burial2(i,j,m,:) & - & + wts2*burial2(i,j,n,:) & - & + wts2*burial(i,j,:) enddo enddo enddo enddo !$OMP END PARALLEL DO + + !$OMP PARALLEL DO PRIVATE(l,i) + do j=1,jj + do l=1,isp(j) + do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) ! time smoothing (analog to tmsmt2.F) + burial2(i,j,m,:) = wts1*burial2(i,j,m,:) & + + wts2*burial2(i,j,n,:) & + + wts2*burial(i,j,:) + enddo + enddo + enddo + !$OMP END PARALLEL DO + !$OMP PARALLEL DO PRIVATE(k,kn,l,i) do k=1,ks kn=k+nns @@ -453,13 +470,22 @@ subroutine hamocc2blom(m,n,mm,nn) do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) sedlay2(i,j,kn,:) = sedlay(i,j,k,:) ! new time level replaces old time level here powtra2(i,j,kn,:) = powtra(i,j,k,:) - burial2(i,j,n,:) = burial(i,j,:) enddo enddo enddo enddo !$OMP END PARALLEL DO + !$OMP PARALLEL DO PRIVATE(l,i) + do j=1,jj + do l=1,isp(j) + do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) + burial2(i,j,n,:) = burial(i,j,:) + enddo + enddo + enddo + !$OMP END PARALLEL DO + endif ! .not. use_sedbypass ! --- ------------------------------------------------------------------ diff --git a/hamocc/mo_ocprod.F90 b/hamocc/mo_ocprod.F90 index 4641e37b..c3fcd2aa 100644 --- a/hamocc/mo_ocprod.F90 +++ b/hamocc/mo_ocprod.F90 @@ -43,11 +43,11 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) ! I.Kriest, *GEOMAR* 2016-08-11 ! - Modified stoichiometry for denitrification (affects NO3, N2, Alk) ! J.Schwinger, *UNI-RESEARCH* 2017-08-30 - ! - Removed split of the layer that only partly falls into the euphotic zone. Loops are + ! - Removed split of the layer that only partly falls into the euphotic zone. Loops are ! now calculated over ! (1) layers that are completely or partly in the euphotoc zone ! (2) layers that do not lie within the euphotic zone. - ! - Moved the accumulation of global fields for output to routine hamocc4bgc. + ! - Moved the accumulation of global fields for output to routine hamocc4bgc. ! The accumulation of local fields has been moved to the end of this routine. ! A.Moree, *GFI, Bergen* 2018-04-12 ! - new version of carbon isotope code @@ -72,7 +72,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) alow1,alow2,alow3,calmax,cellmass, & cellsink,dustd1,dustd2,dustd3,dustsink,fractdim, & fse,fsh,nmldmin,plower,pupper,sinkexp,stick,tmfac, & - tsfac,vsmall,zdis,wmin,wmax,wlin,rbro,bifr13,bifr14, & + tsfac,vsmall,zdis,wmin,wmax,wlin,rbro, & dmsp1,dmsp2,dmsp3,dmsp4,dmsp5,dmsp6,dms_gamma, & fbro1,fbro2,atten_f,atten_c,atten_uv,atten_w,bkopal,bkphy,bkzoo use mo_biomod, only: bsiflx0100,bsiflx0500,bsiflx1000,bsiflx2000,bsiflx4000,bsiflx_bot, & @@ -80,8 +80,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) carflx0100,carflx0500,carflx1000,carflx2000,carflx4000,carflx_bot, & expoor,exposi,expoca,intdnit,intdms_bac,intdmsprod,intdms_uv, & intphosy,int_chbr3_prod,int_chbr3_uv, & - phosy3d,abs_oce,strahl,asize3d,wmass,wnumb,eps3d,bifr13_perm, & - growth_co2 + phosy3d,abs_oce,strahl,asize3d,wmass,wnumb,eps3d use mo_param1_bgc, only: ialkali,ian2o,iano3,icalc,idet,idms,idoc,ifdust, & igasnit,iiron,iopal,ioxygen,iphosph,iphy,isco212, & isilica,izoo,iadust,inos,ibromo, & @@ -124,13 +123,15 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) real :: absorption,absorption_uv real :: dmsprod,dms_bac,dms_uv,dms_ph real :: dtr,dz - real :: wpocd,wcald,wopald,dagg + real :: wpocd,wcald,wopald,wdustd,dagg real :: wcal,wdust,wopal,wpoc ! sedbypass real :: florca,flcaca,flsil ! cisonew real :: phygrowth real :: phosy13,phosy14 + real :: growth_co2 + real :: bifr13,bifr14,bifr13_perm real :: grazing13,grazing14 real :: graton13,graton14 real :: gratpoc13,gratpoc14 @@ -271,7 +272,8 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) !$OMP ,graton13,graton14,gratpoc13,gratpoc14,grawa13,grawa14 & !$OMP ,phosy13,phosy14,bacfra13,bacfra14,phymor13,phymor14,zoomor13 & !$OMP ,zoomor14,excdoc13,excdoc14,exud13,exud14,export13,export14 & - !$OMP ,delcar13,delcar14,dtr13,dtr14,bifr13,bifr14 & + !$OMP ,delcar13,delcar14,dtr13,dtr14,bifr13,bifr14,bifr13_perm & + !$OMP ,growth_co2,phygrowth & !$OMP ,bro_beta,bro_uv & !$OMP ,i,k) @@ -979,7 +981,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) ! C(k,T+dt)=(ddpo(k)*C(k,T)+w*dt*C(k-1,T+dt))/(ddpo(k)+w*dt) ! sedimentation=w*dt*C(ks,T+dt) ! - !$OMP PARALLEL DO PRIVATE(kdonor,wpoc,wpocd,wcal,wcald,wopal,wopald,wnos,wnosd,dagg,i,k) + !$OMP PARALLEL DO PRIVATE(kdonor,wpoc,wpocd,wcal,wcald,wopal,wopald,wdust,wdustd,tco,tcn,q,wnos,wnosd,dagg,i,k) ORDERED do j = 1,kpje do i = 1,kpie @@ -990,7 +992,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) kdonor = 1 do k = 1,kpke - + !$OMP ORDERED ! Sum up total column inventory before sinking scheme if( pddpo(i,j,k) > dp_min ) then tco( 1) = tco( 1) + ocetra(i,j,k,idet )*pddpo(i,j,k) @@ -1025,6 +1027,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) wnos = wnumb(i,j,k) wnosd = wnumb(i,j,kdonor) wdust = dustsink + wdustd = dustsink dagg = dustagg(i,j,k) else if (use_WLIN) then wpoc = min(wmin+wlin*ptiestu(i,j,k), wmax) @@ -1034,6 +1037,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) wopal = wopal_const wopald = wopal_const wdust = wdust_const + wdustd = wdust_const dagg = 0.0 else wpoc = wpoc_const @@ -1043,6 +1047,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) wopal = wopal_const wopald = wopal_const wdust = wdust_const + wdustd = wdust_const dagg = 0.0 endif @@ -1050,6 +1055,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) wpocd = 0.0 wcald = 0.0 wopald = 0.0 + wdustd = 0.0 if (use_AGG) then wnosd = 0.0 else if (use_WLIN) then @@ -1057,15 +1063,15 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) endif endif - ocetra(i,j,k,iopal) = (ocetra(i,j,k, iopal) *pddpo(i,j,k) & - & + ocetra(i,j,kdonor,iopal) *wopald)/ (pddpo(i,j,k)+wopal) - ocetra(i,j,k,ifdust) = (ocetra(i,j,k, ifdust)*pddpo(i,j,k) & - & + ocetra(i,j,kdonor,ifdust)*wdust) / (pddpo(i,j,k)+wdust) & - - dagg ocetra(i,j,k,idet) = (ocetra(i,j,k, idet) *pddpo(i,j,k) & + ocetra(i,j,kdonor,idet) *wpocd) / (pddpo(i,j,k)+wpoc) ocetra(i,j,k,icalc) = (ocetra(i,j,k, icalc) *pddpo(i,j,k) & & + ocetra(i,j,kdonor,icalc) *wcald) / (pddpo(i,j,k)+wcal) + ocetra(i,j,k,iopal) = (ocetra(i,j,k, iopal) *pddpo(i,j,k) & + & + ocetra(i,j,kdonor,iopal) *wopald)/ (pddpo(i,j,k)+wopal) + ocetra(i,j,k,ifdust) = (ocetra(i,j,k, ifdust)*pddpo(i,j,k) & + & + ocetra(i,j,kdonor,ifdust)*wdustd)/ (pddpo(i,j,k)+wdust) & + - dagg if (use_cisonew) then ocetra(i,j,k,idet13) = (ocetra(i,j,k, idet13) *pddpo(i,j,k) & & + ocetra(i,j,kdonor,idet13) *wpocd) / (pddpo(i,j,k)+wpoc) @@ -1137,7 +1143,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) tcn(12) = tcn(12) + ocetra(i,j,k,icalc14)*pddpo(i,j,k) endif endif - + !$OMP END ORDERED enddo ! loop k=1,kpke @@ -1349,7 +1355,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) ! over the water column. Detritus is kept as detritus, while opal and CaCO3 ! are remineralised instantanously - !$OMP PARALLEL DO PRIVATE(dz,florca,flcaca,flsil,flor13,flor14,flca13,flca14,i,k) + !$OMP PARALLEL DO PRIVATE(dz,florca,flcaca,flsil,flor13,flor14,flca13,flca14,i,k) ORDERED do j=1,kpje do i = 1,kpie if(omask(i,j) > 0.5) then @@ -1357,9 +1363,9 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) ! calculate depth of water column dz = 0.0 do k = 1,kpke - + !$OMP ORDERED if( pddpo(i,j,k) > dp_min ) dz = dz+pddpo(i,j,k) - + !$OMP END ORDERED enddo florca = prorca(i,j)/dz @@ -1397,7 +1403,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) endif ! omask > 0.5 enddo enddo - + !$OMP END PARALLEL DO endif ! use_sedbypass if (use_PBGC_OCNP_TIMESTEP) then diff --git a/hamocc/mo_param_bgc.F90 b/hamocc/mo_param_bgc.F90 index 642640b8..3009df95 100644 --- a/hamocc/mo_param_bgc.F90 +++ b/hamocc/mo_param_bgc.F90 @@ -65,7 +65,7 @@ module mo_param_bgc public :: re14to,prei13,prei14,ctochl public :: atten_w,atten_c,atten_uv,atten_f public :: perc_diron,fesoly,phytomi,pi_alpha - public :: dyphy,tf2,tf1,tf0,tff,bifr13,bifr14,c14_t_half + public :: dyphy,tf2,tf1,tf0,tff,bifr13_ini,bifr14_ini,c14_t_half public :: rbro,fbro1,fbro2,grami public :: calmax,remido public :: dustd1,dustd2,dustd3,dustsink @@ -168,8 +168,8 @@ module mo_param_bgc real, protected :: dyphy = 0.004 ! 1/d -mortality rate of phytoplankton ! Initial fractionation during photosynthesis - real :: bifr13 = 0.98 - real :: bifr14 + real, protected :: bifr13_ini = 0.98 + real, protected :: bifr14_ini ! N2-Fixation following the parameterization in Kriest and Oschlies, 2015. ! Factors tf2, tf1 and tf0 are a polynomial (2nd order) @@ -415,7 +415,7 @@ subroutine calc_param_biol() ! AFTER reading the namelist: ! calulate parameters that depend on other tunable parameters ! - bifr14 = bifr13**2 + bifr14_ini = bifr13_ini**2 perc_diron = fetune * 0.035 * 0.01 / 55.85 @@ -591,8 +591,8 @@ subroutine write_parambgc() write(io_stdo_bgc,*) '* atm_c13 = ',atm_c13 write(io_stdo_bgc,*) '* d13C_atm = ',d13C_atm write(io_stdo_bgc,*) '* atm_c14 = ',atm_c14 - write(io_stdo_bgc,*) '* bifr13 = ',bifr13 - write(io_stdo_bgc,*) '* bifr14 = ',bifr14 + write(io_stdo_bgc,*) '* bifr13_ini = ',bifr13_ini + write(io_stdo_bgc,*) '* bifr14_ini = ',bifr14_ini write(io_stdo_bgc,*) '* c14fac = ',c14fac write(io_stdo_bgc,*) '* prei13 = ',prei13 write(io_stdo_bgc,*) '* prei14 = ',prei14 diff --git a/hamocc/mo_powach.F90 b/hamocc/mo_powach.F90 index 80934f5f..f77be13a 100644 --- a/hamocc/mo_powach.F90 +++ b/hamocc/mo_powach.F90 @@ -38,7 +38,8 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) isssc14,issso13,issso14,safediv use mo_carbch, only: co3,keqb,ocetra,sedfluxo use mo_chemcon, only: calcon - use mo_param_bgc, only: rnit,ro2ut,disso_sil,silsat,disso_poc,sed_denit,disso_caco3 + use mo_param_bgc, only: rnit,rcar,rdnit1,rdnit2,ro2ut,disso_sil,silsat,disso_poc,sed_denit, & + disso_caco3 use mo_sedmnt, only: porwat,porsol,powtra,produs,prcaca,prorca,seddw,sedhpl,sedlay, & silpro,pror13,pror14,prca13,prca14 use mo_vgrid, only: kbo,bolay @@ -60,9 +61,9 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) integer :: i,j,k,l real :: sedb1(kpie,0:ks),sediso(kpie,0:ks) real :: solrat(kpie,ks),powcar(kpie,ks) - real :: aerob(kpie,ks),anaerob(kpie,ks) - real :: aerob13(kpie,ks),anaerob13(kpie,ks) ! cisonew - real :: aerob14(kpie,ks),anaerob14(kpie,ks) ! cisonew + real :: aerob(kpie,ks),anaerob(kpie,ks),sulf(kpie,ks) + real :: aerob13(kpie,ks),anaerob13(kpie,ks),sulf13(kpie,ks) ! cisonew + real :: aerob14(kpie,ks),anaerob14(kpie,ks),sulf14(kpie,ks) ! cisonew real :: dissot, undsa, posol real :: umfa, denit, saln, rrho, alk, c, sit, pt real :: K1, K2, Kb, Kw, Ks1, Kf, Ksi, K1p, K2p, K3p @@ -96,11 +97,14 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) powcar(i,k) = 0. anaerob(i,k)= 0. aerob(i,k) = 0. + sulf(i,k) = 0. if (use_cisonew) then anaerob13(i,k)=0. aerob13(i,k) =0. + sulf13(i,k) =0. anaerob14(i,k)=0. aerob14(i,k) =0. + sulf14(i,k) =0. endif enddo enddo @@ -286,7 +290,7 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) do i = 1, kpie if(omask(i,j) > 0.5) then if(powtra(i,j,k,ipowaox) < 1.e-6) then - posol = denit * min(0.5*powtra(i,j,k,ipowno3)/114., sedlay(i,j,k,issso12)) + posol = denit * min(0.25*powtra(i,j,k,ipowno3)/rdnit2, sedlay(i,j,k,issso12)) umfa = porsol(i,j,k)/porwat(i,j,k) anaerob(i,k) = posol*umfa !this has P units: kmol P/m3 of pore water if (use_cisonew) then @@ -299,8 +303,8 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) endif sedlay(i,j,k,issso12) = sedlay(i,j,k,issso12) - posol powtra(i,j,k,ipowaph) = powtra(i,j,k,ipowaph) + posol*umfa - powtra(i,j,k,ipowno3) = powtra(i,j,k,ipowno3) - 98.*posol*umfa - powtra(i,j,k,ipown2) = powtra(i,j,k,ipown2) + 57.*posol*umfa + powtra(i,j,k,ipowno3) = powtra(i,j,k,ipowno3) - rdnit1*posol*umfa + powtra(i,j,k,ipown2) = powtra(i,j,k,ipown2) + rdnit2*posol*umfa if (use_cisonew) then sedlay(i,j,k,issso13) = sedlay(i,j,k,issso13) - poso13 sedlay(i,j,k,issso14) = sedlay(i,j,k,issso14) - poso14 @@ -317,15 +321,14 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) if(powtra(i,j,k,ipowaox) < 3.e-6 .and. powtra(i,j,k,ipowno3) < 3.e-6) then posol = denit * sedlay(i,j,k,issso12) ! remineralization of poc umfa = porsol(i,j,k) / porwat(i,j,k) - !this overwrites anaerob from denitrification. added =anaerob+..., works - anaerob(i,k) = anaerob(i,k) + posol*umfa !this has P units: kmol P/m3 of pore water + sulf(i,k) = posol*umfa !this has P units: kmol P/m3 of pore water if (use_cisonew) then rato13 = sedlay(i,j,k,issso13) / (sedlay(i,j,k,issso12)+safediv) rato14 = sedlay(i,j,k,issso14) / (sedlay(i,j,k,issso12)+safediv) poso13 = posol * rato13 poso14 = posol * rato14 - anaerob13(i,k) = anaerob13(i,k) + poso13*umfa !this has P units: kmol P/m3 of pore water - anaerob14(i,k) = anaerob13(i,k) + poso14*umfa !this has P units: kmol P/m3 of pore water + sulf13(i,k) = poso13*umfa !this has P units: kmol P/m3 of pore water + sulf14(i,k) = poso14*umfa !this has P units: kmol P/m3 of pore water endif sedlay(i,j,k,issso12) = sedlay(i,j,k,issso12) - posol powtra(i,j,k,ipowaph) = powtra(i,j,k,ipowaph) + posol*umfa @@ -352,8 +355,8 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) if(omask(i,j) > 0.5) then saln= min( 40., max( 0., psao(i,j,kbo(i,j)))) rrho= prho(i,j,kbo(i,j)) - alk = (powtra(i,j,k,ipowaal) - (anaerob(i,k)+aerob(i,k))*16.) / rrho - c = (powtra(i,j,k,ipowaic) + (anaerob(i,k)+aerob(i,k))*122.) / rrho + alk = (powtra(i,j,k,ipowaal) - (sulf(i,k)+aerob(i,k))*(rnit+1.) + anaerob(i,k)*(rdnit1-1.)) / rrho + c = (powtra(i,j,k,ipowaic) + (anaerob(i,k)+aerob(i,k)+sulf(i,k))*rcar) / rrho sit = powtra(i,j,k,ipowasi) / rrho pt = powtra(i,j,k,ipowaph) / rrho ah1 = sedhpl(i,j,k) @@ -429,7 +432,7 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) do i = 1, kpie if(omask(i,j) > 0.5) then sedlay(i,j,1,isssc12) = & - & sedlay(i,j,1,isssc12) + prcaca(i,j) / (porsol(i,j,1)*seddw(1)) + & sedlay(i,j,1,isssc12) + prcaca(i,j) / (porsol(i,j,1)*seddw(1)) if (use_cisonew) then sedlay(i,j,1,isssc13) = & & sedlay(i,j,1,isssc13) + prca13(i,j) / (porsol(i,j,1)*seddw(1)) @@ -459,16 +462,16 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,lspin) endif sedlay(i,j,k,isssc12) = sedlay(i,j,k,isssc12) - posol powtra(i,j,k,ipowaic) = powtra(i,j,k,ipowaic) & - + posol * umfa + (aerob(i,k) + anaerob(i,k)) * 122. + & + posol * umfa + (aerob(i,k) + anaerob(i,k) + sulf(i,k)) * rcar powtra(i,j,k,ipowaal) = powtra(i,j,k,ipowaal) & - + 2. * posol * umfa - 16. * (aerob(i,k) + anaerob(i,k)) + & + 2. * posol * umfa - (rnit+1.)*(aerob(i,k) + sulf(i,k)) + (rdnit1-1.)*anaerob(i,k) if (use_cisonew) then sedlay(i,j,k,isssc13) = sedlay(i,j,k,isssc13) - poso13 sedlay(i,j,k,isssc14) = sedlay(i,j,k,isssc14) - poso14 powtra(i,j,k,ipowc13) = powtra(i,j,k,ipowc13) + poso13 * umfa & - + (aerob13(i,k) + anaerob13(i,k)) * 122. + & + (aerob13(i,k) + anaerob13(i,k) + sulf13(i,k)) * rcar powtra(i,j,k,ipowc14) = powtra(i,j,k,ipowc14) + poso14 * umfa & - + (aerob14(i,k) + anaerob14(i,k)) * 122. + & + (aerob14(i,k) + anaerob14(i,k) + sulf14(i,k)) * rcar endif endif enddo diff --git a/hamocc/mo_vgrid.F90 b/hamocc/mo_vgrid.F90 index e395cd49..a585444d 100644 --- a/hamocc/mo_vgrid.F90 +++ b/hamocc/mo_vgrid.F90 @@ -90,8 +90,9 @@ subroutine set_vgrid(kpie,kpje,kpke,pddpo) ! --- depth of layer kpke+1 centre ptiestu(:,:,kpke+1)=9000. - !$OMP PARALLEL DO PRIVATE(j,i) + do k=1,kpke + !$OMP PARALLEL DO PRIVATE(j,i) do j=1,kpje do i=1,kpie @@ -102,8 +103,9 @@ subroutine set_vgrid(kpie,kpje,kpke,pddpo) enddo enddo + !$OMP END PARALLEL DO enddo - !$OMP END PARALLEL DO + kbo(:,:) =1 bolay(:,:)=0.0