Merge branch 'master' into LGM

f90/cfo.f90

@@ -244,7 +244,7 @@ CONTAINS
     EI = zoo%E*zoo%IC
     zoo%ngr = (EI*(1E0_dp - zoo%ca) - zoo%poc*(zoo%cf*zoo%Af + Rm))*fQ
     zoo%RC  = EI - zoo%ngr           ! respiration
-    zoo%XC = zoo%IC*(1E0_dp - zoo%E) ! egestion
+    zoo%XC = zoo%IC - EI             ! egestion
   END SUBROUTINE zxrex
 
   ! egest extra C in food

f90/cmo.f90

@@ -336,14 +336,22 @@ CONTAINS
     CLASS(ocm), INTENT(INOUT) :: phy
     REAL(dp), INTENT(INOUT) :: RC
     REAL(dp), INTENT(IN) :: phin, zN
-    REAL(dp) :: paravg, PARb, aim, depfrac
+    REAL(dp) :: paravg, PARb, aim, depfrac, l1, l2
     IF ((phy%PARs.GT.phy%I0).AND.(phy%mu0.GT.1E-30_dp)) THEN            ! (A4)
        PARb = MAX(phy%I0, phy%PARb)
        depfrac = LOG(phy%PARs/PARb)
        paravg = (phy%PARs - PARb)/depfrac
        aim = phy%alpha*paravg/phy%mu0
-       phy%tch = 1E0_dp/phy%zC + (1E0_dp - lambertw((1E0_dp + RC/(phy%daylen*phy%mu0))&
-            *EXP(1E0_dp + MIN(aim/phy%zC, 2E2_dp)), 0, 0))/aim
+       IF (aim.GT.7.e2_dp*phy%zC) THEN
+          ! approximation of lambertw(x) for large x: W0(x) = l1 - l2 + l2/l1
+          ! with l1 = ln(x), l2 = ln(l1)
+          l1 = LOG(1.E0_dp + RC/(phy%daylen*phy%mu0)) + 1E0_dp + aim/phy%zC
+          l2 = LOG(l1)
+          phy%tch = 1.E0_dp/phy%zC + (1E0_dp - l1 + l2 - l2/l1)/aim
+       ELSE
+          phy%tch = 1.E0_dp/phy%zC + (1E0_dp - lambertw((1.E0_dp + RC/(phy%daylen*phy%mu0))&
+               *EXP(1E0_dp + aim/phy%zC), 0, 0))/aim
+       END IF
        phy%SI = physi(phy, phin, phy%PARs, PARb, depfrac)
        phy%theta = phy%tch*phy%fC
        RC = RC*depfrac/phin          ! for calculation of phy%A in SR grow

source/mom/tracer.F

@@ -390,7 +390,7 @@
      &           *1.2
 #  endif
             dayfrac = min( 1., -tan(tlat(i,jrow)/radian)*tan(declin))
-            dayfrac = max(1e-12, acos(max(-1., dayfrac))/pi)
+            dayfrac = max(1.e-12, acos(max(-1., dayfrac))/pi)
 #  if defined O_embm
 ! swr is in W/m2 = kg/s3, dnswr is in g/s3 = mW/m2
 ! from Andreas: swr is 24h-avg.
@@ -399,7 +399,7 @@
             swr = 200.
 #  endif
 #  if defined O_opem
-     &            /dayfrac ! with cmo, swr is mean DAYTIME irradiance
+     &            /MAX(dayfrac, 1.e-3) ! with cmo, swr is mean DAYTIME irradiance
             glb = swr
             expo_c = 0.0
             expo_p = 0.0
@@ -781,23 +781,23 @@
                 sg_bdeni = (0.06 + 0.19*0.99
      &                        **(max(t(i,k,j,io2,taum1)*1000,trcmin)
      &                        - max(t(i,k,j,ino3,taum1),trcmin)))
-#   if defined O_opem
+#     if defined O_opem
      &               *max(expo_c*sgb,trcmin)
-#   else
+#     else
      &               *max(expo*sgb,trcmin)*6.625
-#   endif
+#     endif
                 sg_bdeni = min(sg_bdeni, sgb*expo)
                 sg_bdeni = max(sg_bdeni, 0.)
 
                 sg_bdeni = sg_bdeni*(0.5 + lno3)*no3flag*din15flag
 
-#   if defined O_LGM
+#     if defined O_LGM
                 if (k .eq. 1) then
                    sg_bdeni = 0.0
                 else if (k .eq. 2) then
                    sg_bdeni = 0.08*sg_bdeni
                 endif
-#   endif
+#     endif
                 snpzd(5) = snpzd(5) + sgb*expo - sg_bdeni
 
 #    else
@@ -846,9 +846,21 @@
                 endif
 #   else
                 if (t(i,k,j,io2,taum1) .lt. o2min) then
+#    if defined O_LGM
+! remineralise only the fraction of Fe as the assumed existence of the current grid cell
+! this is 0 for k.eq.1 and 0.08 for k.eq.2
+                   if (k .eq. 2) then
+                      snpzd(itdfe) = snpzd(itdfe) + expofe*sgb*0.08
+                      bfe = bfe + expofe*sgb*0.08
+                      expofe = expofe - sgb*expofe*0.08
+                   elseif (k .gt. 2) then
+#    endif
                    snpzd(itdfe) = snpzd(itdfe) + expofe*sgb
                    bfe = bfe + expofe*sgb
                    expofe = expofe - sgb*expofe
+#    if defined O_LGM
+                   endif
+#    endif
                 endif
 #   endif
 
@@ -905,11 +917,11 @@
                 bdeni = (0.06 + 0.19*0.99
      &               **(max(t(i,k,j,io2,taum1)*1000,trcmin)
      &               - max(t(i,k,j,ino3,taum1),trcmin)))
-#   if defined O_opem
+#     if defined O_opem
      &               *max(expo_c,trcmin)
-#   else
+#     else
      &               *max(expo,trcmin)*6.625
-#   endif
+#     endif
                 bdeni = min(bdeni, expo)
                 bdeni = max(bdeni, 0.)
 
@@ -917,6 +929,13 @@ c                if (k .ge. 1 .and. k .le. 10)
 c     &                  bdeni = sgbdfac*bdeni
                 bdeni = bdeni*(0.5 + lno3)*no3flag*din15flag
 
+#      if defined O_LGM
+                if (k .eq. 1) then
+                   bdeni = 0.0
+                else if (k .eq. 2) then
+                   bdeni = 0.08*bdeni
+                endif
+#      endif
                 snpzd(5) = snpzd(5) + expo - bdeni
 #    else
                 snpzd(5) = snpzd(5) + expo
@@ -927,13 +946,31 @@ c     &                  bdeni = sgbdfac*bdeni
                 fesed = fetopsed*bct*expo_p
 #    else
                 fesed = fetopsed*bct*redptn*expo
+#    endif
+! As for bdeni, also reduce sedimentary iron release from the sediment.
+#    if defined O_LGM
+                if (k .eq. 1) then
+                   fesed = 0.0
+                else if (k .eq. 2) then
+                   fesed = fesed*0.08
+                endif
 #    endif
                 bfe = fesed
                 snpzd(itdfe) = snpzd(itdfe) + fesed
                 if (t(i,k,j,io2,taum1) .lt. o2min) then
+#    if defined O_LGM
+                   if (k .eq. 2) then
+                      snpzd(itdfe) = snpzd(itdfe) + expofe*0.08
+                      bfe = bfe + expofe*0.08
+                      expofe = expofe - expofe*0.08
+                   elseif (k .gt. 2) then
+#    endif
                    snpzd(itdfe) = snpzd(itdfe) + expofe
                    bfe = bfe + expofe
                    expofe = 0.
+#    if defined O_LGM
+                   endif
+#    endif
                 endif
                 rremife(i,k,j) = rremife(i,k,j) + bfe
 #    if defined O_npzd_iron_diagnostics