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.
seddta.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/SED/seddta.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 seddta
2   !!======================================================================
3   !!                     ***  MODULE  seddta  ***
4   !! Sediment data  :  read sediment input data from a file
5   !!=====================================================================
6
7   !! * Modules used
8   USE sed
9   USE sedarr
10   USE phycst, ONLY : rday
11   USE iom
12   USE lib_mpp         ! distribued memory computing library
13
14   IMPLICIT NONE
15   PRIVATE
16
17   !! * Routine accessibility
18   PUBLIC sed_dta   !
19
20   !! *  Module variables
21   REAL(wp) ::  rsecday  ! number of second per a day
22   REAL(wp) ::  conv2    ! [kg/m2/month]-->[g/cm2/s] ( 1 month has 30 days )
23
24   !! * Substitutions
25#  include "do_loop_substitute.h90"
26   !! $Id$
27CONTAINS
28
29   !!---------------------------------------------------------------------------
30   !!   sed_dta  : read the NetCDF data file in online version using module iom
31   !!---------------------------------------------------------------------------
32
33   SUBROUTINE sed_dta( kt, Kbb, Kmm )
34      !!----------------------------------------------------------------------
35      !!                   ***  ROUTINE sed_dta  ***
36      !!                   
37      !! ** Purpose :   Reads data from a netcdf file and
38      !!                initialization of rain and pore water (k=1) components
39      !!
40      !!
41      !!   History :
42      !!        !  04-10  (N. Emprin, M. Gehlen )  Original code
43      !!        !  06-04  (C. Ethe)  Re-organization ; Use of iom
44      !!----------------------------------------------------------------------
45
46      !! Arguments
47      INTEGER, INTENT(in) ::  kt         ! time-step
48      INTEGER, INTENT(in) ::  Kbb, Kmm   ! time level indices
49
50      !! * Local declarations
51      INTEGER  ::  ji, jj, js, jw, ikt
52
53      REAL(wp), DIMENSION(jpoce) :: zdtap, zdtag
54      REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3
55      REAL(wp) :: zf0, zf1, zf2, zkapp, zratio, zdep
56
57      !----------------------------------------------------------------------
58
59      ! Initialization of sediment variable
60      ! Spatial dimension is merged, and unity converted if needed
61      !-------------------------------------------------------------
62
63      IF( ln_timing )  CALL timing_start('sed_dta')
64
65      IF (lwp) THEN
66         WRITE(numsed,*)
67         WRITE(numsed,*) ' sed_dta : Bottom layer fields'
68         WRITE(numsed,*) ' ~~~~~~'
69         WRITE(numsed,*) ' Data from SMS model'
70         WRITE(numsed,*)
71      ENDIF
72
73
74      ! open file
75      IF( kt == nitsed000 ) THEN
76         IF (lwp) WRITE(numsed,*) ' sed_dta : Sediment fields'
77         dtsed = r2dttrc
78         rsecday = 60.* 60. * 24.
79!         conv2   = 1.0e+3 / ( 1.0e+4 * rsecday * 30. )
80         conv2 = 1.0e+3 /  1.0e+4 
81         rdtsed(2:jpksed) = dtsed / ( denssol * por1(2:jpksed) )
82      ENDIF
83
84      ! Initialization of temporaries arrays 
85      zdtap(:)    = 0. 
86      zdtag(:)    = 0. 
87
88      ! reading variables
89      IF (lwp) WRITE(numsed,*)
90      IF (lwp) WRITE(numsed,*) ' sed_dta : Bottom layer fields at time  kt = ', kt
91      ! reading variables
92      !
93      !    Sinking speeds of detritus is increased with depth as shown
94      !    by data and from the coagulation theory
95      !    -----------------------------------------------------------
96      IF (ln_sediment_offline) THEN
97         DO_2D_11_11
98            ikt = mbkt(ji,jj)
99            zwsbio4(ji,jj) = wsbio2 / rday
100            zwsbio3(ji,jj) = wsbio  / rday
101         END_2D
102      ELSE
103         DO_2D_11_11
104            ikt = mbkt(ji,jj)
105            zdep = e3t(ji,jj,ikt,Kmm) / r2dttrc
106            zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) / rday )
107            zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) / rday )
108         END_2D
109      ENDIF
110
111      trc_data(:,:,:) = 0.
112      DO_2D_11_11
113         ikt = mbkt(ji,jj)
114         IF ( tmask(ji,jj,ikt) == 1 ) THEN
115            trc_data(ji,jj,1)   = tr(ji,jj,ikt,jpsil,Kbb)
116            trc_data(ji,jj,2)   = tr(ji,jj,ikt,jpoxy,Kbb)
117            trc_data(ji,jj,3)   = tr(ji,jj,ikt,jpdic,Kbb)
118            trc_data(ji,jj,4)   = tr(ji,jj,ikt,jpno3,Kbb) / 7.625
119            trc_data(ji,jj,5)   = tr(ji,jj,ikt,jppo4,Kbb) / 122.
120            trc_data(ji,jj,6)   = tr(ji,jj,ikt,jptal,Kbb)
121            trc_data(ji,jj,7)   = tr(ji,jj,ikt,jpnh4,Kbb) / 7.625
122            trc_data(ji,jj,8)   = 0.0
123            trc_data(ji,jj,9)   = 28.0E-3
124            trc_data(ji,jj,10)  = tr(ji,jj,ikt,jpfer,Kbb)
125            trc_data(ji,jj,11 ) = MIN(tr(ji,jj,ikt,jpgsi,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3
126            trc_data(ji,jj,12 ) = MIN(tr(ji,jj,ikt,jppoc,Kbb), 1E-4) * zwsbio3(ji,jj) * 1E3
127            trc_data(ji,jj,13 ) = MIN(tr(ji,jj,ikt,jpgoc,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3
128            trc_data(ji,jj,14)  = MIN(tr(ji,jj,ikt,jpcal,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3
129            trc_data(ji,jj,15)  = ts(ji,jj,ikt,jp_tem,Kmm)
130            trc_data(ji,jj,16)  = ts(ji,jj,ikt,jp_sal,Kmm)
131            trc_data(ji,jj,17 ) = ( tr(ji,jj,ikt,jpsfe,Kbb) * zwsbio3(ji,jj) + tr(ji,jj,ikt,jpbfe,Kbb)  &
132            &                     * zwsbio4(ji,jj)  ) * 1E3 / ( trc_data(ji,jj,12 ) + trc_data(ji,jj,13 ) + rtrn )
133            trc_data(ji,jj,17 ) = MIN(1E-3, trc_data(ji,jj,17 ) )
134         ENDIF
135      END_2D
136
137      ! Pore water initial concentration [mol/l] in  k=1
138      !-------------------------------------------------
139      DO jw = 1, jpwat
140         CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jw), trc_data(1:jpi,1:jpj,jw), iarroce(1:jpoce) )
141      END DO
142      !  Solid components :
143      !-----------------------
144      !  Sinking fluxes for OPAL in mol.m-2.s-1 ; conversion in mol.cm-2.s-1
145      CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsopal), trc_data(1:jpi,1:jpj,11), iarroce(1:jpoce) ) 
146      rainrm_dta(1:jpoce,jsopal) = rainrm_dta(1:jpoce,jsopal) * 1e-4
147      !  Sinking fluxes for POC in mol.m-2.s-1 ; conversion in mol.cm-2.s-1
148      CALL pack_arr ( jpoce, zdtap(1:jpoce), trc_data(1:jpi,1:jpj,12) , iarroce(1:jpoce) )     
149      CALL pack_arr ( jpoce, zdtag(1:jpoce), trc_data(1:jpi,1:jpj,13) , iarroce(1:jpoce) )
150      DO ji = 1, jpoce
151!        zkapp  = MIN( (1.0 - 0.02 ) * reac_poc, 3731.0 * max(100.0, zkbot(ji) )**(-1.011) / ( 365.0 * 24.0 * 3600.0 ) )
152!        zkapp   = MIN( 0.98 * reac_poc, 100.0 * max(100.0, zkbot(ji) )**(-0.6) / ( 365.0 * 24.0 * 3600.0 ) )
153!        zratio = ( ( 1.0 - 0.02 ) * reac_poc + 0.02 * reac_poc * 0. - zkapp) / ( ( 0.02 - 1.0 ) * reac_poc / 100. - 0.02 * reac_poc * 0. + zkapp )
154!        zf1    = ( 0.02 * (reac_poc - reac_poc * 0.) + zkapp - reac_poc ) / ( reac_poc / 100. - reac_poc )
155!        zf1    = MIN(0.98, MAX(0., zf1 ) )
156         zf1    = 0.48
157         zf0    = 1.0 - 0.02 - zf1
158         zf2    = 0.02
159         rainrm_dta(ji,jspoc) =   ( zdtap(ji) +  zdtag(ji) ) * 1e-4 * zf0
160         rainrm_dta(ji,jspos) =   ( zdtap(ji) +  zdtag(ji) ) * 1e-4 * zf1
161         rainrm_dta(ji,jspor) =   ( zdtap(ji) +  zdtag(ji) ) * 1e-4 * zf2
162      END DO
163      !  Sinking fluxes for Calcite in mol.m-2.s-1 ; conversion in mol.cm-2.s-1
164      CALL pack_arr ( jpoce,  rainrm_dta(1:jpoce,jscal), trc_data(1:jpi,1:jpj,14), iarroce(1:jpoce) )
165      rainrm_dta(1:jpoce,jscal) = rainrm_dta(1:jpoce,jscal) * 1e-4
166      ! vector temperature [°C] and salinity
167      CALL pack_arr ( jpoce,  temp(1:jpoce), trc_data(1:jpi,1:jpj,15), iarroce(1:jpoce) )
168      CALL pack_arr ( jpoce,  salt(1:jpoce), trc_data(1:jpi,1:jpj,16), iarroce(1:jpoce) )
169     
170      ! Clay rain rate in [mol/(cm**2.s)]
171      ! inputs data in [kg.m-2.sec-1] ---> 1e+3/(1e+4) [g.cm-2.s-1]   
172      ! divided after by molecular weight g.mol-1     
173      CALL pack_arr ( jpoce,  rainrm_dta(1:jpoce,jsclay), dust(1:jpi,1:jpj), iarroce(1:jpoce) )
174      rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * conv2 / mol_wgt(jsclay)   &
175      &                            + wacc(1:jpoce) * por1(2) * denssol / mol_wgt(jsclay) / ( rsecday * 365.0 )
176      rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * 0.965
177      rainrm_dta(1:jpoce,jsfeo)  = rainrm_dta(1:jpoce,jsclay) * mol_wgt(jsclay) / mol_wgt(jsfeo) * 0.035 / 0.965
178!    rainrm_dta(1:jpoce,jsclay) = 1.0E-4 * conv2 / mol_wgt(jsclay)
179
180      ! Iron monosulphide rain rates. Set to 0
181      rainrm_dta(1:jpoce,jsfes)  = 0. 
182
183      ! Fe/C ratio in sinking particles that fall to the sediments
184      CALL pack_arr ( jpoce,  fecratio(1:jpoce), trc_data(1:jpi,1:jpj,17), iarroce(1:jpoce) )
185
186      sedligand(:,1) = 1.E-9
187
188      ! sediment pore water at 1st layer (k=1)
189      DO jw = 1, jpwat
190         pwcp(1:jpoce,1,jw) = pwcp_dta(1:jpoce,jw)
191      ENDDO
192
193      !  rain
194      DO js = 1, jpsol
195         rainrm(1:jpoce,js) = rainrm_dta(1:jpoce,js)
196      ENDDO
197
198      ! Calculation of raintg of each sol. comp.: rainrm in [g/(cm**2.s)]
199      DO js = 1, jpsol
200         rainrg(1:jpoce,js) = rainrm(1:jpoce,js) *  mol_wgt(js)
201      ENDDO
202
203      ! Calculation of raintg = total massic flux rained in each cell (sum of sol. comp.)
204      raintg(:) = 0.
205      DO js = 1, jpsol
206         raintg(1:jpoce) = raintg(1:jpoce) + rainrg(1:jpoce,js)
207      ENDDO
208
209      ! computation of dzdep = total thickness of solid material rained [cm] in each cell
210      dzdep(1:jpoce) = raintg(1:jpoce) * rdtsed(2) 
211
212      IF( lk_iomput ) THEN
213          IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", dust(:,:) * conv2 * 1E4 )
214          IF( iom_use("sflxcal" ) )  CALL iom_put( "sflxcal", trc_data(:,:,13) )
215          IF( iom_use("sflxbsi" ) )  CALL iom_put( "sflxbsi", trc_data(:,:,10) )
216          IF( iom_use("sflxpoc" ) )  CALL iom_put( "sflxpoc", trc_data(:,:,11) + trc_data(:,:,12) )
217      ENDIF
218
219      IF( ln_timing )  CALL timing_stop('sed_dta')
220     
221   END SUBROUTINE sed_dta
222
223END MODULE seddta
Note: See TracBrowser for help on using the repository browser.