diff --git a/hamocc/beleg_vars.F90 b/hamocc/beleg_vars.F90 index 22cd8dc4..f7d68963 100644 --- a/hamocc/beleg_vars.F90 +++ b/hamocc/beleg_vars.F90 @@ -220,18 +220,16 @@ SUBROUTINE BELEG_VARS(kpaufr,kpie,kpje,kpke,kbnd,pddpo,prho,omask, & ! Initialise preformed tracers in the mixed layer; note that the ! whole field has been initialised to zero above - DO k=1,kmle DO j=1,kpje DO i=1,kpie IF(omask(i,j) .GT. 0.5) THEN - ocetra(i,j,k,iprefo2) =ocetra(i,j,k,ioxygen) - ocetra(i,j,k,iprefpo4)=ocetra(i,j,k,iphosph) - ocetra(i,j,k,iprefalk)=ocetra(i,j,k,ialkali) - ocetra(i,j,k,iprefdic)=ocetra(i,j,k,isco212) + ocetra(i,j,1:kmle(i,j),iprefo2) = ocetra(i,j,1:kmle(i,j),ioxygen) + ocetra(i,j,1:kmle(i,j),iprefpo4) = ocetra(i,j,1:kmle(i,j),iphosph) + ocetra(i,j,1:kmle(i,j),iprefalk) = ocetra(i,j,1:kmle(i,j),ialkali) + ocetra(i,j,1:kmle(i,j),iprefdic) = ocetra(i,j,1:kmle(i,j),isco212) ENDIF ENDDO ENDDO - ENDDO ! Initial values for sediment diff --git a/hamocc/carchm.F90 b/hamocc/carchm.F90 index d86013b8..bab04daf 100644 --- a/hamocc/carchm.F90 +++ b/hamocc/carchm.F90 @@ -96,11 +96,11 @@ SUBROUTINE CARCHM(kpie,kpje,kpke,kbnd, & !********************************************************************** use mo_carbch, only: atm,atmflx,co2fxd,co2fxu,co2star,co3,hi,keqb,kwco2sol,ocetra,omegaa,omegac,pco2d,satn2o,satoxy use mo_chemcon, only: al1,al2,al3,al4,an0,an1,an2,an3,an4,an5,an6,atn2o,bl1,bl2,bl3,calcon,ox0,ox1,ox2,ox3,ox4,ox5,ox6, & - & oxyco,tzero - use mo_control_bgc, only: dtbgc + & oxyco,tzero + use mo_control_bgc, only: dtbgc use mo_param1_bgc, only: ialkali,iatmo2,iatmco2,iatmdms,iatmn2,iatmn2o,ian2o,icalc,idicsat,idms,igasnit,ioxygen,iphosph, & - & isco212,isilica - use mo_vgrid, only: dp_min,kbo,ptiestu + & isco212,isilica + use mo_vgrid, only: dp_min,kmle,kbo,ptiestu #ifdef BROMO use mo_param1_bgc, only: iatmbromo,ibromo @@ -390,8 +390,7 @@ SUBROUTINE CARCHM(kpie,kpje,kpke,kbnd, & ta = ocetra(i,j,k,ialkali) / rrho CALL carchm_solve_DICsat(s,atco2*rpp0,ta,sit,pt,Kh,K1,K2,Kb,Kw,Ks1,Kf, & Ksi,K1p,K2p,K3p,tc_sat,niter) - ocetra(i,j,k, idicsat)=tc_sat * rrho ! convert mol/kg to kmol/m^3 - ocetra(i,j,k+1,idicsat)=tc_sat * rrho ! k+1 = the rest of the mixed layer + ocetra(i,j,1:kmle(i,j),idicsat) = tc_sat * rrho ! convert mol/kg to kmlo/m^3 #ifdef cisonew ! Ocean-Atmosphere fluxes for carbon isotopes diff --git a/hamocc/cyano.F90 b/hamocc/cyano.F90 index e39fa678..f3f696df 100644 --- a/hamocc/cyano.F90 +++ b/hamocc/cyano.F90 @@ -17,7 +17,7 @@ ! along with BLOM. If not, see https://www.gnu.org/licenses/. - SUBROUTINE CYANO(kpie,kpje,kpke,kbnd,pddpo,omask,ptho) +SUBROUTINE CYANO(kpie,kpje,kpke,kbnd,pddpo,omask,ptho) !********************************************************************** ! !**** *CYANO* - . @@ -61,74 +61,69 @@ SUBROUTINE CYANO(kpie,kpje,kpke,kbnd,pddpo,omask,ptho) ! . !********************************************************************** - use mo_carbch, only: ocetra - use mo_biomod, only: bluefix,intnfix,rnit,tf0,tf1,tf2,tff - use mo_param1_bgc, only: ialkali,iano3,igasnit,iphosph,ioxygen - use mo_vgrid, only: kmle + use mo_carbch, only: ocetra + use mo_biomod, only: bluefix,intnfix,rnit,tf0,tf1,tf2,tff + use mo_param1_bgc, only: ialkali,iano3,igasnit,iphosph,ioxygen + use mo_vgrid, only: kmle #ifdef natDIC - use mo_param1_bgc, only: inatalkali + use mo_param1_bgc, only: inatalkali #endif - implicit none + implicit none - INTEGER, intent(in) :: kpie,kpje,kpke,kbnd - REAL, intent(in) :: pddpo(kpie,kpje,kpke) - REAL, intent(in) :: omask(kpie,kpje) - REAL, intent(in) :: ptho(1-kbnd:kpie+kbnd,1-kbnd:kpje+kbnd,kpke) + INTEGER, intent(in) :: kpie,kpje,kpke,kbnd + REAL, intent(in) :: pddpo(kpie,kpje,kpke) + REAL, intent(in) :: omask(kpie,kpje) + REAL, intent(in) :: ptho(1-kbnd:kpie+kbnd,1-kbnd:kpje+kbnd,kpke) - ! Local variables - INTEGER :: i,j,k - REAL :: oldocetra,dano3 - REAL :: ttemp,nfixtfac + ! Local variables + INTEGER :: i,j,k + REAL :: oldocetra,dano3 + REAL :: ttemp,nfixtfac + + intnfix(:,:)=0.0 - intnfix(:,:)=0.0 - ! -! N-fixation by cyano bacteria (followed by remineralisation and nitrification), +! N-fixation by cyano bacteria (followed by remineralisation and nitrification), ! it is assumed here that this process is limited to the mixed layer ! - DO k=1,kmle -!$OMP PARALLEL DO PRIVATE(i,oldocetra,dano3,ttemp,nfixtfac) - DO j=1,kpje - DO i=1,kpie - IF(omask(i,j).gt.0.5) THEN - IF(ocetra(i,j,k,iano3).LT.(rnit*ocetra(i,j,k,iphosph))) THEN + DO j=1,kpje + DO i=1,kpie + IF(omask(i,j).gt.0.5) THEN + DO k=1,kmle(i,j) + IF(ocetra(i,j,k,iano3).LT.(rnit*ocetra(i,j,k,iphosph))) THEN - oldocetra = ocetra(i,j,k,iano3) - ttemp = min(40.,max(-3.,ptho(i,j,k))) + oldocetra = ocetra(i,j,k,iano3) + ttemp = min(40.,max(-3.,ptho(i,j,k))) ! Temperature dependence of nitrogen fixation, Kriest and Oschlies 2015. - nfixtfac = MAX(0.0,tf2*ttemp*ttemp + tf1*ttemp + tf0)/tff + nfixtfac = MAX(0.0,tf2*ttemp*ttemp + tf1*ttemp + tf0)/tff - ocetra(i,j,k,iano3)=ocetra(i,j,k,iano3)*(1-bluefix*nfixtfac) & - & +bluefix*nfixtfac*rnit*ocetra(i,j,k,iphosph) + ocetra(i,j,k,iano3)=ocetra(i,j,k,iano3)*(1-bluefix*nfixtfac) & + & + bluefix*nfixtfac*rnit*ocetra(i,j,k,iphosph) - dano3=ocetra(i,j,k,iano3)-oldocetra + dano3=ocetra(i,j,k,iano3)-oldocetra - ocetra(i,j,k,igasnit)=ocetra(i,j,k,igasnit)-dano3*(1./2.) + ocetra(i,j,k,igasnit)=ocetra(i,j,k,igasnit)-dano3*(1./2.) ! Note: to fix one mole N2 requires: N2+H2O+y*O2 = 2* HNO3 <-> y=2.5 mole O2. ! I.e., to release one mole HNO3 = H+ + NO3- requires 1.25 mole O2 - ocetra(i,j,k,ioxygen)=ocetra(i,j,k,ioxygen)-dano3*1.25 + ocetra(i,j,k,ioxygen)=ocetra(i,j,k,ioxygen)-dano3*1.25 ! Nitrogen fixation followed by remineralisation and nitrification decreases ! alkalinity by 1 mole per mole nitrogen fixed (Wolf-Gladrow et al. 2007) - ocetra(i,j,k,ialkali)=ocetra(i,j,k,ialkali)-dano3 + ocetra(i,j,k,ialkali)=ocetra(i,j,k,ialkali)-dano3 #ifdef natDIC - ocetra(i,j,k,inatalkali)=ocetra(i,j,k,inatalkali)-dano3 + ocetra(i,j,k,inatalkali)=ocetra(i,j,k,inatalkali)-dano3 #endif - intnfix(i,j) = intnfix(i,j) + & - & (ocetra(i,j,k,iano3)-oldocetra)*pddpo(i,j,k) - - ENDIF - ENDIF - ENDDO - ENDDO -!$OMP END PARALLEL DO - ENDDO - + intnfix(i,j) = intnfix(i,j) + & + & (ocetra(i,j,k,iano3)-oldocetra)*pddpo(i,j,k) + ENDIF + ENDDO + ENDIF + ENDDO + ENDDO - RETURN - END +END SUBROUTINE CYANO diff --git a/hamocc/mo_apply_rivin.F90 b/hamocc/mo_apply_rivin.F90 index 697990d1..0dfbc528 100644 --- a/hamocc/mo_apply_rivin.F90 +++ b/hamocc/mo_apply_rivin.F90 @@ -120,28 +120,28 @@ subroutine apply_rivin(kpie,kpje,kpke,pddpo,omask,rivin) ! Distribute riverine inputs over the model mixed layer volij = 0. - DO k=1,kmle + DO k=1,kmle(i,j) volij=volij+pddpo(i,j,k) ENDDO ! DIC is updated using the assumtions that a_t=a_c+a_n and DIC=a_c (a_t: total ! alkalinity, a_c: carbonate alkalinity, a_n: contribution of nutrients to a_t). - ocetra(i,j,1:kmle,iano3) = ocetra(i,j,1:kmle,iano3) + rivin(i,j,irdin)*fdt/volij - ocetra(i,j,1:kmle,iphosph) = ocetra(i,j,1:kmle,iphosph) + rivin(i,j,irdip)*fdt/volij - ocetra(i,j,1:kmle,isilica) = ocetra(i,j,1:kmle,isilica) + rivin(i,j,irsi) *fdt/volij - ocetra(i,j,1:kmle,isco212) = ocetra(i,j,1:kmle,isco212) + rivin(i,j,iralk)*fdt/volij & + ocetra(i,j,1:kmle(i,j),iano3) = ocetra(i,j,1:kmle(i,j),iano3) + rivin(i,j,irdin)*fdt/volij + ocetra(i,j,1:kmle(i,j),iphosph) = ocetra(i,j,1:kmle(i,j),iphosph) + rivin(i,j,irdip)*fdt/volij + ocetra(i,j,1:kmle(i,j),isilica) = ocetra(i,j,1:kmle(i,j),isilica) + rivin(i,j,irsi) *fdt/volij + ocetra(i,j,1:kmle(i,j),isco212) = ocetra(i,j,1:kmle(i,j),isco212) + rivin(i,j,iralk)*fdt/volij & + rivin(i,j,irdin)*fdt/volij & + rivin(i,j,irdip)*fdt/volij - ocetra(i,j,1:kmle,ialkali) = ocetra(i,j,1:kmle,ialkali) + rivin(i,j,iralk)*fdt/volij + ocetra(i,j,1:kmle(i,j),ialkali) = ocetra(i,j,1:kmle(i,j),ialkali) + rivin(i,j,iralk)*fdt/volij #ifdef natDIC - ocetra(i,j,1:kmle,inatsco212) = ocetra(i,j,1:kmle,inatsco212) + rivin(i,j,iralk)*fdt/volij & + ocetra(i,j,1:kmle(i,j),inatsco212) = ocetra(i,j,1:kmle(i,j),inatsco212) + rivin(i,j,iralk)*fdt/volij & + rivin(i,j,irdin)*fdt/volij & + rivin(i,j,irdip)*fdt/volij - ocetra(i,j,1:kmle,inatalkali) = ocetra(i,j,1:kmle,inatalkali) + rivin(i,j,iralk)*fdt/volij + ocetra(i,j,1:kmle(i,j),inatalkali) = ocetra(i,j,1:kmle(i,j),inatalkali) + rivin(i,j,iralk)*fdt/volij #endif - ocetra(i,j,1:kmle,iiron) = ocetra(i,j,1:kmle,iiron) + rivin(i,j,iriron)*fdt/volij*dFe_frac - ocetra(i,j,1:kmle,idoc) = ocetra(i,j,1:kmle,idoc) + rivin(i,j,irdoc)*fdt/volij - ocetra(i,j,1:kmle,idet) = ocetra(i,j,1:kmle,idet) + rivin(i,j,irdet)*fdt/volij + ocetra(i,j,1:kmle(i,j),iiron) = ocetra(i,j,1:kmle(i,j),iiron) + rivin(i,j,iriron)*fdt/volij*dFe_frac + ocetra(i,j,1:kmle(i,j),idoc) = ocetra(i,j,1:kmle(i,j),idoc) + rivin(i,j,irdoc)*fdt/volij + ocetra(i,j,1:kmle(i,j),idet) = ocetra(i,j,1:kmle(i,j),idet) + rivin(i,j,irdet)*fdt/volij rivinflx(i,j,irdin) = rivin(i,j,irdin)*fdt rivinflx(i,j,irdip) = rivin(i,j,irdip)*fdt diff --git a/hamocc/mo_intfcblom.F90 b/hamocc/mo_intfcblom.F90 index 68227f7b..e0d78b3b 100644 --- a/hamocc/mo_intfcblom.F90 +++ b/hamocc/mo_intfcblom.F90 @@ -244,14 +244,16 @@ subroutine blom2hamocc(m,n,mm,nn) !****************************************************************************** ! use mod_constants, only: onem - use mod_xc, only: ii,jdm,jj,kdm,kk,ifp,isp,ilp,idm + use mod_xc, only: ii,jdm,jj,kdm,kk,ifp,isp,ilp,idm use mod_grid, only: scpx,scpy use mod_state, only: dp,temp,saln use mod_eos, only: rho,p_alpha + use mod_difest, only: hOBL use mod_tracers, only: ntrbgc,itrbgc,trc use mo_param1_bgc, only: ks,nsedtra,npowtra,natm use mo_carbch, only: ocetra,atm use mo_sedmnt, only: sedlay,powtra,sedhpl,burial + use mo_vgrid, only: kmle, kmle_static implicit none @@ -292,6 +294,11 @@ subroutine blom2hamocc(m,n,mm,nn) ! --- - dimension of grid box in meters bgc_dx(i,j) = scpx(i,j)/1.e2 bgc_dy(i,j) = scpy(i,j)/1.e2 +! +! --- - index of level above OBL depth +! --- isopycninc coords: hOBL(i,j) = hOBL_static = 3. => kmle(i,j) = 2 +! --- hybrid coords: hOBL defined according to cvmix_kpp_compute_kOBL_depth + kmle(i,j) = nint(hOBL(i,j))-1 enddo enddo !$OMP END PARALLEL DO diff --git a/hamocc/mo_vgrid.F90 b/hamocc/mo_vgrid.F90 index 0f7cc08b..e010e92e 100644 --- a/hamocc/mo_vgrid.F90 +++ b/hamocc/mo_vgrid.F90 @@ -53,16 +53,18 @@ module mo_vgrid !****************************************************************************** implicit none - INTEGER, PARAMETER :: kmle = 2 ! k-end index for layers that - ! represent the mixed layer in BLOM + INTEGER, PARAMETER :: kmle_static = 2 ! k-end index for layers that + ! represent the mixed layer in BLOM. + ! Default value used for isopycnic coordinates. REAL, PARAMETER :: dp_ez = 100.0 ! depth of euphotic zone - REAL, PARAMETER :: dp_min = 1.0E-12 ! min layer thickness layers thinner + REAL, PARAMETER :: dp_min = 1.0E-12 ! min layer thickness layers thinner ! than this are ignored by HAMOCC REAL, PARAMETER :: dp_min_sink = 1.0 ! min layer thickness for sinking (layers thinner than ! this are ignored and set to the concentration of the ! layer above). Note that the bottom layer index kbo(i,j) ! is defined as the lowermost layer thicker than dp_min_sink. + INTEGER, DIMENSION(:,:), ALLOCATABLE :: kmle INTEGER, DIMENSION(:,:), ALLOCATABLE :: kbo INTEGER, DIMENSION(:,:), ALLOCATABLE :: kwrbioz INTEGER, DIMENSION(:,:), ALLOCATABLE :: k0100,k0500,k1000,k2000,k4000 @@ -263,6 +265,17 @@ subroutine alloc_mem_vgrid(kpie,kpje,kpke) ptiestw(:,:,:) = 0.0 + IF(mnproc.eq.1) THEN + WRITE(io_stdo_bgc,*)'Memory allocation for variable kmle ...' + WRITE(io_stdo_bgc,*)'First dimension : ',kpie + WRITE(io_stdo_bgc,*)'Second dimension : ',kpje + ENDIF + + ALLOCATE(kmle(kpie,kpje),stat=errstat) + if(errstat.ne.0) stop 'not enough memory kmle' + kmle(:,:) = kmle_static + + IF(mnproc.eq.1) THEN WRITE(io_stdo_bgc,*)'Memory allocation for variable kbo ...' WRITE(io_stdo_bgc,*)'First dimension : ',kpie diff --git a/hamocc/ocprod.F90 b/hamocc/ocprod.F90 index 2677b483..d98bc43c 100644 --- a/hamocc/ocprod.F90 +++ b/hamocc/ocprod.F90 @@ -938,7 +938,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) ! Set minimum particle number to nmldmin in the mixed layer. This is to prevent ! very small values of nos (and asscociated high sinking speed if there is mass) ! in high latitudes during winter - if ( k <= kmle ) then + if ( k <= kmle(i,j) ) then ocetra(i,j,k,inos) = MAX(nmldmin,ocetra(i,j,k,inos)) endif @@ -974,7 +974,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph) ! As a first step, assume that shear in the mixed layer is high and ! zero below. - if ( k <= kmle ) then + if ( k <= kmle(i,j) ) then fshear = fsh else fshear = 0. diff --git a/hamocc/preftrc.F90 b/hamocc/preftrc.F90 index 34a9161c..a33280d1 100644 --- a/hamocc/preftrc.F90 +++ b/hamocc/preftrc.F90 @@ -16,7 +16,7 @@ ! along with BLOM. If not, see https://www.gnu.org/licenses/. - SUBROUTINE PREFTRC(kpie,kpje,omask) +SUBROUTINE PREFTRC(kpie,kpje,omask) !**************************************************************** ! !**** *PREFTRC* - update preformed tracers in the mixed layer. @@ -43,31 +43,27 @@ SUBROUTINE PREFTRC(kpie,kpje,omask) ! !************************************************************************** - use mo_carbch, only: ocetra - use mo_param1_bgc, only: ialkali,ioxygen,iphosph,iprefalk,iprefdic,iprefo2,iprefpo4,isco212 - use mo_vgrid, only: kmle + use mo_carbch, only: ocetra + use mo_param1_bgc, only: ialkali,ioxygen,iphosph,iprefalk,iprefdic,iprefo2,iprefpo4,isco212 + use mo_vgrid, only: kmle - implicit none + implicit none - INTEGER :: kpie,kpje - REAL :: omask(kpie,kpje) + INTEGER :: kpie,kpje + REAL :: omask(kpie,kpje) - INTEGER :: i,j,k + INTEGER :: i,j - do k=1,kmle -!$OMP PARALLEL DO PRIVATE(i) - do j=1,kpje - do i=1,kpie + do j=1,kpje + do i=1,kpie if (omask(i,j) .gt. 0.5 ) then - ocetra(i,j,k,iprefo2) =ocetra(i,j,k,ioxygen) - ocetra(i,j,k,iprefpo4)=ocetra(i,j,k,iphosph) - ocetra(i,j,k,iprefalk)=ocetra(i,j,k,ialkali) - ocetra(i,j,k,iprefdic)=ocetra(i,j,k,isco212) + ocetra(i,j,1:kmle(i,j),iprefo2) = ocetra(i,j,1:kmle(i,j),ioxygen) + ocetra(i,j,1:kmle(i,j),iprefpo4) = ocetra(i,j,1:kmle(i,j),iphosph) + ocetra(i,j,1:kmle(i,j),iprefalk) = ocetra(i,j,1:kmle(i,j),ialkali) + ocetra(i,j,1:kmle(i,j),iprefdic) = ocetra(i,j,1:kmle(i,j),isco212) endif - enddo - enddo -!$OMP END PARALLEL DO - enddo + enddo + enddo - END SUBROUTINE PREFTRC +END SUBROUTINE PREFTRC diff --git a/phy/mod_difest.F b/phy/mod_difest.F index 2d2dde96..5bd4b2b9 100644 --- a/phy/mod_difest.F +++ b/phy/mod_difest.F @@ -78,6 +78,10 @@ module mod_difest implicit none c private +c +c Initialize hOBL with hOBL_static = 3. for consistency with bulk +c mixed layer formulation in iHAMOCC: kmle = nint(hOBL) - 1 = 2 + real, PARAMETER :: hOBL_static = 3. c real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm+1) :: . rig @@ -85,6 +89,8 @@ module mod_difest . du2l,drhol,up,vp real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: . OBLdepth + real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: + . hOBL integer, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm) :: . mskv,msku integer, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: @@ -174,7 +180,7 @@ module mod_difest . cs=98.96,minOBLdepth=1.0) c public :: OBLdepth, inivar_difest, init_difest, difest_isobml, - . difest_lateral_hybrid, difest_vertical_hybrid + . difest_lateral_hybrid, difest_vertical_hybrid, hOBL c contains c @@ -203,6 +209,7 @@ subroutine inivar_difest enddo do i=1-nbdy,ii+nbdy OBLdepth(i,j)=spval + hOBL(i,j) = hOBL_static enddo enddo c$OMP END PARALLEL DO @@ -904,7 +911,6 @@ subroutine difest_vertical_hyb(m,n,mm,nn,k1m,k1n) real :: Simmons_coeff, zBottomMinusOffset real :: bl1, bl2, bl3, bl4 integer ki, ksfc, ktmp, kOBL, kn1 - real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: hOBL c surf_layer_ext = 0.1 bl1 = 8e-5