Fix potential OMP issues and bug fix in wrapping the burial fields (#273

) * fix missing OMP END * fix potential race condition * Fix multiple time smoothing for burial fields and ommit multiple rewriting of burial fields * Prevent from potential race condition, and make water column fields private * ORDERED clause for sedbypass=True case
NorESMhub · Sep 26, 2023 · 68845e1 · 68845e1
1 parent 8690509
commit 68845e1
Show file tree

Hide file tree

Showing 4 changed files with 48 additions and 19 deletions.
diff --git a/hamocc/dipowa.F90 b/hamocc/dipowa.F90
@@ -84,8 +84,7 @@ subroutine dipowa(kpie,kpje,kpke,omask,lspin)
 !ik needed for boundary layer ventilation in fast sediment routine
   real :: bolven(kpie)                 ! bottom layer ventilation rate
 
-!$OMP PARALLEL DO                            &
-!$OMP&PRIVATE(i,k,iv,l,bolven,tredsy,sedb1,aprior,iv_oc)
+!$OMP PARALLEL DO PRIVATE(i,k,iv,l,bolven,tredsy,sedb1,aprior,iv_oc)
   j_loop: do j=1,kpje
 
 ! calculate bottom ventilation rate for scaling of sediment-water exchange
@@ -206,5 +205,5 @@ 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_intfcblom.F90 b/hamocc/mo_intfcblom.F90
@@ -377,12 +377,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
 
 ! --- ------------------------------------------------------------------
@@ -486,15 +496,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
@@ -503,12 +523,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
 
 ! --- ------------------------------------------------------------------

diff --git a/hamocc/mo_vgrid.F90 b/hamocc/mo_vgrid.F90
@@ -110,8 +110,8 @@ 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
 
@@ -122,8 +122,8 @@ subroutine set_vgrid(kpie,kpje,kpke,pddpo)
 
     enddo
     enddo
-  enddo
 !$OMP END PARALLEL DO
+  enddo
 
 
   kbo(:,:)  =1

diff --git a/hamocc/ocprod.F90 b/hamocc/ocprod.F90
@@ -131,7 +131,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph)
   real,    intent(in) :: ptho(1-kbnd:kpie+kbnd,1-kbnd:kpje+kbnd,kpke)
   real,    intent(in) :: pi_ph(kpie,kpje)
 
-  ! Local varaibles
+  ! Local variables
   integer :: i,j,k,l
   integer :: is,kdonor
   integer, parameter :: nsinkmax = 12
@@ -1049,11 +1049,11 @@ 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,wdust,wdustd   &
+!$OMP PARALLEL DO PRIVATE(kdonor,wpoc,wpocd,wcal,wcald,wopal,wopald,wdust,wdustd,tco,tcn,q    &
 #if defined(AGG)
-!$OMP ,wnos,wnosd,dagg                                                &
+!$OMP ,wnos,wnosd,dagg                                                                        &
 #endif
-!$OMP ,i,k)
+!$OMP ,i,k) ORDERED
   do j = 1,kpje
   do i = 1,kpie
 
@@ -1064,7 +1064,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)
@@ -1217,7 +1217,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
 
 
@@ -1439,17 +1439,17 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph)
 #ifdef cisonew
 !$OMP ,flor13,flor14,flca13,flca14                                      &
 #endif
-!$OMP ,i,k)
+!$OMP ,i,k) ORDERED
   do j=1,kpje
   do i = 1,kpie
      if(omask(i,j) > 0.5) then
 
         ! 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