New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
p4zsink.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zsink.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.
  • Property svn:keywords set to Id
File size: 9.4 KB
Line 
1MODULE p4zsink
2   !!======================================================================
3   !!                         ***  MODULE p4zsink  ***
4   !! TOP :  PISCES  vertical flux of particulate matter due to gravitational sinking
5   !!======================================================================
6   !! History :   1.0  !  2004     (O. Aumont) Original code
7   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
8   !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Change aggregation formula
9   !!             3.5  !  2012-07  (O. Aumont) Introduce potential time-splitting
10   !!----------------------------------------------------------------------
11   !!   p4z_sink       :  Compute vertical flux of particulate matter due to gravitational sinking
12   !!   p4z_sink_init  :  Unitialisation of sinking speed parameters
13   !!   p4z_sink_alloc :  Allocate sinking speed variables
14   !!----------------------------------------------------------------------
15   USE oce_trc         !  shared variables between ocean and passive tracers
16   USE trc             !  passive tracers common variables
17   USE sms_pisces      !  PISCES Source Minus Sink variables
18   USE trcsink         !  General routine to compute sedimentation
19   USE prtctl_trc      !  print control for debugging
20   USE iom             !  I/O manager
21   USE lib_mpp
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   p4z_sink         ! called in p4zbio.F90
27   PUBLIC   p4z_sink_init    ! called in trcsms_pisces.F90
28   PUBLIC   p4z_sink_alloc
29
30   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinking, sinking2  !: POC sinking fluxes
31   !                                                          !  (different meanings depending on the parameterization)
32   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkingn, sinking2n  !: POC sinking fluxes
33   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkingp, sinking2p  !: POC sinking fluxes
34   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkcal, sinksil   !: CaCO3 and BSi sinking fluxes
35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkfer            !: Small BFe sinking fluxes
36   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkfer2           !: Big iron sinking fluxes
37
38   INTEGER  :: ik100
39
40   !! * Substitutions
41#  include "do_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
44   !! $Id$
45   !! Software governed by the CeCILL license (see ./LICENSE)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   !!----------------------------------------------------------------------
50   !!   'standard sinking parameterisation'                  ???
51   !!----------------------------------------------------------------------
52
53   SUBROUTINE p4z_sink ( kt, knt, Kbb, Kmm, Krhs )
54      !!---------------------------------------------------------------------
55      !!                     ***  ROUTINE p4z_sink  ***
56      !!
57      !! ** Purpose :   Compute vertical flux of particulate matter due to
58      !!                gravitational sinking
59      !!
60      !! ** Method  : - ???
61      !!---------------------------------------------------------------------
62      INTEGER, INTENT(in) :: kt, knt
63      INTEGER, INTENT(in) :: Kbb, Kmm, Krhs  ! time level indices
64      INTEGER  ::   ji, jj, jk
65      CHARACTER (len=25) :: charout
66      REAL(wp) :: zmax, zfact
67      !!---------------------------------------------------------------------
68      !
69      IF( ln_timing )   CALL timing_start('p4z_sink')
70
71      ! Initialization of some global variables
72      ! ---------------------------------------
73      prodpoc(:,:,:) = 0.
74      conspoc(:,:,:) = 0.
75      prodgoc(:,:,:) = 0.
76      consgoc(:,:,:) = 0.
77
78      !
79      !    Sinking speeds of detritus is increased with depth as shown
80      !    by data and from the coagulation theory
81      !    -----------------------------------------------------------
82      DO_3D_11_11( 1, jpkm1 )
83         zmax  = MAX( heup_01(ji,jj), hmld(ji,jj) )
84         zfact = MAX( 0., gdepw(ji,jj,jk+1,Kmm) - zmax ) / wsbio2scale
85         wsbio4(ji,jj,jk) = wsbio2 + MAX(0., ( wsbio2max - wsbio2 )) * zfact
86      END_3D
87
88      ! limit the values of the sinking speeds to avoid numerical instabilities 
89      wsbio3(:,:,:) = wsbio
90
91      !
92      !  Initializa to zero all the sinking arrays
93      !   -----------------------------------------
94      sinking (:,:,:) = 0.e0
95      sinking2(:,:,:) = 0.e0
96      sinkcal (:,:,:) = 0.e0
97      sinkfer (:,:,:) = 0.e0
98      sinksil (:,:,:) = 0.e0
99      sinkfer2(:,:,:) = 0.e0
100
101      !   Compute the sedimentation term using p4zsink2 for all the sinking particles
102      !   -----------------------------------------------------
103      CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinking , jppoc, rfact2 )
104      CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinkfer , jpsfe, rfact2 )
105      CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinking2, jpgoc, rfact2 )
106      CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinkfer2, jpbfe, rfact2 )
107      CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinksil , jpgsi, rfact2 )
108      CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinkcal , jpcal, rfact2 )
109
110      IF( ln_p5z ) THEN
111         sinkingn (:,:,:) = 0.e0
112         sinking2n(:,:,:) = 0.e0
113         sinkingp (:,:,:) = 0.e0
114         sinking2p(:,:,:) = 0.e0
115
116         !   Compute the sedimentation term using p4zsink2 for all the sinking particles
117         !   -----------------------------------------------------
118         CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinkingn , jppon, rfact2 )
119         CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinkingp , jppop, rfact2 )
120         CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinking2n, jpgon, rfact2 )
121         CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinking2p, jpgop, rfact2 )
122      ENDIF
123
124     ! Total carbon export per year
125     IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
126        &   t_oce_co2_exp = glob_sum( 'p4zsink', ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * e1e2t(:,:) * tmask(:,:,1) )
127     !
128     IF( lk_iomput .AND.  knt == nrdttrc ) THEN
129       zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
130       !
131       CALL iom_put( "EPC100"  , ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * zfact * tmask(:,:,1) ) ! Export of carbon at 100m
132       CALL iom_put( "EPFE100" , ( sinkfer(:,:,ik100) + sinkfer2(:,:,ik100) ) * zfact * tmask(:,:,1) ) ! Export of iron at 100m
133       CALL iom_put( "EPCAL100", sinkcal(:,:,ik100) * zfact * tmask(:,:,1) )      ! Export of calcite at 100m
134       CALL iom_put( "EPSI100" , sinksil(:,:,ik100) * zfact * tmask(:,:,1) )          ! Export of bigenic silica at 100m
135       CALL iom_put( "EXPC"    , ( sinking(:,:,:) + sinking2(:,:,:) ) * zfact * tmask(:,:,:) ) ! Export of carbon in the water column
136       CALL iom_put( "EXPFE"   , ( sinkfer(:,:,:) + sinkfer2(:,:,:) ) * zfact * tmask(:,:,:) ) ! Export of iron 
137       CALL iom_put( "EXPCAL"  , sinkcal(:,:,:) * zfact * tmask(:,:,:) )      ! Export of calcite
138       CALL iom_put( "EXPSI"   , sinksil(:,:,:) * zfact * tmask(:,:,:) )      ! Export of bigenic silica
139       CALL iom_put( "tcexp"   , t_oce_co2_exp * zfact )   ! molC/s
140       !
141      ENDIF
142      !
143      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
144         WRITE(charout, FMT="('sink')")
145         CALL prt_ctl_trc_info(charout)
146         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
147      ENDIF
148      !
149      IF( ln_timing )   CALL timing_stop('p4z_sink')
150      !
151   END SUBROUTINE p4z_sink
152
153
154   SUBROUTINE p4z_sink_init
155      !!----------------------------------------------------------------------
156      !!                  ***  ROUTINE p4z_sink_init  ***
157      !!----------------------------------------------------------------------
158      INTEGER :: jk
159      !!----------------------------------------------------------------------
160      !
161      ik100 = 10        !  last level where depth less than 100 m
162      DO jk = jpkm1, 1, -1
163         IF( gdept_1d(jk) > 100. )  ik100 = jk - 1
164      END DO
165      IF (lwp) WRITE(numout,*)
166      IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ',  ik100 + 1
167      IF (lwp) WRITE(numout,*)
168      !
169      t_oce_co2_exp = 0._wp
170      !
171   END SUBROUTINE p4z_sink_init
172
173   INTEGER FUNCTION p4z_sink_alloc()
174      !!----------------------------------------------------------------------
175      !!                     ***  ROUTINE p4z_sink_alloc  ***
176      !!----------------------------------------------------------------------
177      INTEGER :: ierr(2)
178      !!----------------------------------------------------------------------
179      !
180      ierr(:) = 0
181      !
182      ALLOCATE( sinking(jpi,jpj,jpk) , sinking2(jpi,jpj,jpk)                    ,     &               
183         &      sinkcal(jpi,jpj,jpk) , sinksil (jpi,jpj,jpk)                    ,     &               
184         &      sinkfer2(jpi,jpj,jpk)                                           ,     &               
185         &      sinkfer(jpi,jpj,jpk)                                            , STAT=ierr(1) )               
186         !
187      IF( ln_p5z    ) ALLOCATE( sinkingn(jpi,jpj,jpk), sinking2n(jpi,jpj,jpk)   ,     &
188         &                      sinkingp(jpi,jpj,jpk), sinking2p(jpi,jpj,jpk)   , STAT=ierr(2) )
189      !
190      p4z_sink_alloc = MAXVAL( ierr )
191      IF( p4z_sink_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_sink_alloc : failed to allocate arrays.' )
192      !
193   END FUNCTION p4z_sink_alloc
194   
195   !!======================================================================
196END MODULE p4zsink
Note: See TracBrowser for help on using the repository browser.