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.
trcsbc.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/TRP – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/TRP/trcsbc.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: 10.0 KB
Line 
1MODULE trcsbc
2   !!==============================================================================
3   !!                       ***  MODULE  trcsbc  ***
4   !! Ocean passive tracers:  surface boundary condition
5   !!======================================================================
6   !! History :  8.2  !  1998-10  (G. Madec, G. Roullet, M. Imbard)  Original code
7   !!            8.2  !  2001-02  (D. Ludicone)  sea ice and free surface
8   !!            8.5  !  2002-06  (G. Madec)  F90: Free form and module
9   !!            9.0  !  2004-03  (C. Ethe)  adapted for passive tracers
10   !!                 !  2006-08  (C. Deltel) Diagnose ML trends for passive tracers
11   !!==============================================================================
12#if defined key_top
13   !!----------------------------------------------------------------------
14   !!   'key_top'                                                TOP models
15   !!----------------------------------------------------------------------
16   !!   trc_sbc      : update the tracer trend at ocean surface
17   !!----------------------------------------------------------------------
18   USE oce_trc         ! ocean dynamics and active tracers variables
19   USE trc             ! ocean  passive tracers variables
20   USE prtctl_trc      ! Print control for debbuging
21   USE iom
22   USE trd_oce
23   USE trdtra
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   trc_sbc   ! routine called by step.F90
29
30   !! * Substitutions
31#  include "vectopt_loop_substitute.h90"
32#  include "do_loop_substitute.h90"
33   !!----------------------------------------------------------------------
34   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
35   !! $Id$
36   !! Software governed by the CeCILL license (see ./LICENSE)
37   !!----------------------------------------------------------------------
38CONTAINS
39
40   SUBROUTINE trc_sbc ( kt, Kmm, ptr, Krhs )
41      !!----------------------------------------------------------------------
42      !!                  ***  ROUTINE trc_sbc  ***
43      !!                   
44      !! ** Purpose :   Compute the tracer surface boundary condition trend of
45      !!      (concentration/dilution effect) and add it to the general
46      !!       trend of tracer equations.
47      !!
48      !! ** Method :
49      !!      * concentration/dilution effect:
50      !!            The surface freshwater flux modify the ocean volume
51      !!         and thus the concentration of a tracer as :
52      !!            tr(Krhs) = tr(Krhs) + emp * tr(Kmm) / e3t   for k=1
53      !!         where emp, the surface freshwater budget (evaporation minus
54      !!         precipitation ) given in kg/m2/s is divided
55      !!         by 1035 kg/m3 (density of ocean water) to obtain m/s.
56      !!
57      !! ** Action  : - Update the 1st level of tr(:,:,:,:,Krhs) with the trend associated
58      !!                with the tracer surface boundary condition
59      !!
60      !!----------------------------------------------------------------------
61      INTEGER,                                    INTENT(in   ) :: kt        ! ocean time-step index
62      INTEGER,                                    INTENT(in   ) :: Kmm, Krhs ! time level indices
63      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr       ! passive tracers and RHS of tracer equation
64      !
65      INTEGER  ::   ji, jj, jn                      ! dummy loop indices
66      REAL(wp) ::   zse3t, zrtrn, zfact     ! local scalars
67      REAL(wp) ::   zftra, zdtra, ztfx, ztra   !   -      -
68      CHARACTER (len=22) :: charout
69      REAL(wp), DIMENSION(jpi,jpj)   ::   zsfx
70      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrtrd
71      !!---------------------------------------------------------------------
72      !
73      IF( ln_timing )   CALL timing_start('trc_sbc')
74      !
75      ! Allocate temporary workspace
76      IF( l_trdtrc )  ALLOCATE( ztrtrd(jpi,jpj,jpk) )
77      !
78      zrtrn = 1.e-15_wp
79
80      IF( kt == nittrc000 ) THEN
81         IF(lwp) WRITE(numout,*)
82         IF(lwp) WRITE(numout,*) 'trc_sbc : Passive tracers surface boundary condition'
83         IF(lwp) WRITE(numout,*) '~~~~~~~ '
84         !
85         IF( ln_rsttr .AND. .NOT.ln_top_euler .AND.   &                     ! Restart: read in restart  file
86            iom_varid( numrtr, 'sbc_'//TRIM(ctrcnm(1))//'_b', ldstop = .FALSE. ) > 0 ) THEN
87            IF(lwp) WRITE(numout,*) '          nittrc000-1 surface tracer content forcing fields read in the restart file'
88            zfact = 0.5_wp
89            DO jn = 1, jptra
90               CALL iom_get( numrtr, jpdom_autoglo, 'sbc_'//TRIM(ctrcnm(jn))//'_b', sbc_trc_b(:,:,jn) )   ! before tracer content sbc
91            END DO
92         ELSE                                         ! No restart or restart not found: Euler forward time stepping
93           zfact = 1._wp
94           sbc_trc_b(:,:,:) = 0._wp
95         ENDIF
96      ELSE                                         ! Swap of forcing fields
97         IF( ln_top_euler ) THEN
98            zfact = 1._wp
99            sbc_trc_b(:,:,:) = 0._wp
100         ELSE
101            zfact = 0.5_wp
102            sbc_trc_b(:,:,:) = sbc_trc(:,:,:)
103         ENDIF
104         !
105      ENDIF
106
107      ! Coupling online : river runoff is added to the horizontal divergence (hdiv) in the subroutine sbc_rnf_div
108      ! one only consider the concentration/dilution effect due to evaporation minus precipitation + freezing/melting of sea-ice
109      ! Coupling offline : runoff are in emp which contains E-P-R
110      !
111      IF( .NOT.ln_linssh ) THEN  ! online coupling with vvl
112         zsfx(:,:) = 0._wp
113      ELSE                                      ! online coupling free surface or offline with free surface
114         zsfx(:,:) = emp(:,:)
115      ENDIF
116
117      ! 0. initialization
118      SELECT CASE ( nn_ice_tr )
119
120      CASE ( -1 ) ! No tracers in sea ice (null concentration in sea ice)
121         !
122         DO jn = 1, jptra
123            DO_2D_01_00
124               sbc_trc(ji,jj,jn) = zsfx(ji,jj) * r1_rau0 * ptr(ji,jj,1,jn,Kmm)
125            END_2D
126         END DO
127         !
128      CASE ( 0 )  ! Same concentration in sea ice and in the ocean
129         !
130         DO jn = 1, jptra
131            DO_2D_01_00
132               sbc_trc(ji,jj,jn) = ( zsfx(ji,jj) + fmmflx(ji,jj) ) * r1_rau0 * ptr(ji,jj,1,jn,Kmm)
133            END_2D
134         END DO
135         !
136      CASE ( 1 )  ! Specific treatment of sea ice fluxes with an imposed concentration in sea ice
137         !
138         DO jn = 1, jptra
139            DO_2D_01_00
140               zse3t = 1. / e3t(ji,jj,1,Kmm)
141               ! tracer flux at the ice/ocean interface (tracer/m2/s)
142               zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice
143               !                                         ! only used in the levitating sea ice case
144               ! tracer flux only       : add concentration dilution term in net tracer flux, no F-M in volume flux
145               ! tracer and mass fluxes : no concentration dilution term in net tracer flux, F-M term in volume flux
146               ztfx  = zftra                        ! net tracer flux
147               !
148               zdtra = r1_rau0 * ( ztfx + ( zsfx(ji,jj) + fmmflx(ji,jj) ) * ptr(ji,jj,1,jn,Kmm) ) 
149               IF ( zdtra < 0. ) THEN
150                  zdtra  = MAX(zdtra, -ptr(ji,jj,1,jn,Kmm) * e3t(ji,jj,1,Kmm) / r2dttrc )   ! avoid negative concentrations to arise
151               ENDIF
152               sbc_trc(ji,jj,jn) =  zdtra 
153            END_2D
154         END DO
155      END SELECT
156      !
157      CALL lbc_lnk( 'trcsbc', sbc_trc(:,:,:), 'T', 1. )
158      !                                       Concentration dilution effect on tracers due to evaporation & precipitation
159      DO jn = 1, jptra
160         !
161         IF( l_trdtrc )   ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs)  ! save trends
162         !
163         DO_2D_01_00
164            zse3t = zfact / e3t(ji,jj,1,Kmm)
165            ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + ( sbc_trc_b(ji,jj,jn) + sbc_trc(ji,jj,jn) ) * zse3t
166         END_2D
167         !
168         IF( l_trdtrc ) THEN
169            ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs) - ztrtrd(:,:,:)
170            CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_nsr, ztrtrd )
171         END IF
172         !                                                       ! ===========
173      END DO                                                     ! tracer loop
174      !                                                          ! ===========
175      !
176      !                                           Write in the tracer restar  file
177      !                                          *******************************
178      IF( lrst_trc .AND. .NOT.ln_top_euler ) THEN
179         IF(lwp) WRITE(numout,*)
180         IF(lwp) WRITE(numout,*) 'sbc : ocean surface tracer content forcing fields written in tracer restart file ',   &
181            &                    'at it= ', kt,' date= ', ndastp
182         IF(lwp) WRITE(numout,*) '~~~~'
183         DO jn = 1, jptra
184            CALL iom_rstput( kt, nitrst, numrtw, 'sbc_'//TRIM(ctrcnm(jn))//'_b', sbc_trc(:,:,jn) )
185         END DO
186      ENDIF
187      !
188      IF( sn_cfctl%l_prttrc )   THEN
189         WRITE(charout, FMT="('sbc ')") ;  CALL prt_ctl_trc_info(charout)
190                                           CALL prt_ctl_trc( tab4d=ptr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
191      ENDIF
192      IF( l_trdtrc )  DEALLOCATE( ztrtrd )
193      !
194      IF( ln_timing )   CALL timing_stop('trc_sbc')
195      !
196   END SUBROUTINE trc_sbc
197
198#else
199   !!----------------------------------------------------------------------
200   !!   Dummy module :                      NO passive tracer
201   !!----------------------------------------------------------------------
202   USE par_oce
203   USE par_trc
204CONTAINS
205   SUBROUTINE trc_sbc ( kt, Kmm, ptr, Krhs )      ! Empty routine
206      INTEGER,                                    INTENT(in   ) :: kt        ! ocean time-step index
207      INTEGER,                                    INTENT(in   ) :: Kmm, Krhs ! time level indices
208      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr       ! passive tracers and RHS of tracer equation
209      WRITE(*,*) 'trc_sbc: You should not have seen this print! error?', kt
210   END SUBROUTINE trc_sbc
211#endif
212   
213   !!======================================================================
214END MODULE trcsbc
Note: See TracBrowser for help on using the repository browser.