Hamocc hybrid coord2 (NorESMhub#179)

Make the surface mixed layer depth fractional index `hOBL` available for use in iHAMOCC, and adjust the internal iHAMOCC index `kmle` according to `hOBL`. Default value `kmle = 2` is retained for consistency with isopycnic coordinates.

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

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

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

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

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

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

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.

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