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.
sbcssr.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/sbcssr.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: 13.5 KB
Line 
1MODULE sbcssr
2   !!======================================================================
3   !!                       ***  MODULE  sbcssr  ***
4   !! Surface module :  heat and fresh water fluxes a restoring term toward observed SST/SSS
5   !!======================================================================
6   !! History :  3.0  !  2006-06  (G. Madec)  Original code
7   !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   sbc_ssr       : add to sbc a restoring term toward SST/SSS climatology
12   !!   sbc_ssr_init  : initialisation of surface restoring
13   !!----------------------------------------------------------------------
14   USE oce            ! ocean dynamics and tracers
15   USE dom_oce        ! ocean space and time domain
16   USE sbc_oce        ! surface boundary condition
17   USE phycst         ! physical constants
18   USE sbcrnf         ! surface boundary condition : runoffs
19   !
20   USE fldread        ! read input fields
21   USE in_out_manager ! I/O manager
22   USE iom            ! I/O manager
23   USE lib_mpp        ! distribued memory computing library
24   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
25   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   sbc_ssr        ! routine called in sbcmod
31   PUBLIC   sbc_ssr_init   ! routine called in sbcmod
32   PUBLIC   sbc_ssr_alloc  ! routine called in sbcmod
33
34   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   erp   !: evaporation damping   [kg/m2/s]
35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qrp   !: heat flux damping        [w/m2]
36   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   coefice   !: under ice relaxation coefficient
37
38   !                                   !!* Namelist namsbc_ssr *
39   INTEGER, PUBLIC ::   nn_sstr         ! SST/SSS restoring indicator
40   INTEGER, PUBLIC ::   nn_sssr         ! SST/SSS restoring indicator
41   REAL(wp)        ::   rn_dqdt         ! restoring factor on SST and SSS
42   REAL(wp)        ::   rn_deds         ! restoring factor on SST and SSS
43   LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term
44   REAL(wp)        ::   rn_sssr_bnd     ! ABS(Max./Min.) value of erp term [mm/day]
45   INTEGER         ::   nn_sssr_ice     ! Control of restoring under ice
46
47   REAL(wp) , ALLOCATABLE, DIMENSION(:) ::   buffer   ! Temporary buffer for exchange
48   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_sst   ! structure of input SST (file informations, fields read)
49   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_sss   ! structure of input SSS (file informations, fields read)
50
51   !! * Substitutions
52#  include "do_loop_substitute.h90"
53   !!----------------------------------------------------------------------
54   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
55   !! $Id$
56   !! Software governed by the CeCILL license (see ./LICENSE)
57   !!----------------------------------------------------------------------
58CONTAINS
59
60   SUBROUTINE sbc_ssr( kt )
61      !!---------------------------------------------------------------------
62      !!                     ***  ROUTINE sbc_ssr  ***
63      !!
64      !! ** Purpose :   Add to heat and/or freshwater fluxes a damping term
65      !!                toward observed SST and/or SSS.
66      !!
67      !! ** Method  : - Read namelist namsbc_ssr
68      !!              - Read observed SST and/or SSS
69      !!              - at each nscb time step
70      !!                   add a retroaction term on qns    (nn_sstr = 1)
71      !!                   add a damping term on sfx        (nn_sssr = 1)
72      !!                   add a damping term on emp        (nn_sssr = 2)
73      !!---------------------------------------------------------------------
74      INTEGER, INTENT(in   ) ::   kt   ! ocean time step
75      !!
76      INTEGER  ::   ji, jj   ! dummy loop indices
77      REAL(wp) ::   zerp     ! local scalar for evaporation damping
78      REAL(wp) ::   zqrp     ! local scalar for heat flux damping
79      REAL(wp) ::   zsrp     ! local scalar for unit conversion of rn_deds factor
80      REAL(wp) ::   zerp_bnd ! local scalar for unit conversion of rn_epr_max factor
81      INTEGER  ::   ierror   ! return error code
82      !!
83      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
84      TYPE(FLD_N) ::   sn_sst, sn_sss        ! informations about the fields to be read
85      !!----------------------------------------------------------------------
86      !
87      IF( nn_sstr + nn_sssr /= 0 ) THEN
88         !
89         IF( nn_sstr == 1)   CALL fld_read( kt, nn_fsbc, sf_sst )   ! Read SST data and provides it at kt
90         IF( nn_sssr >= 1)   CALL fld_read( kt, nn_fsbc, sf_sss )   ! Read SSS data and provides it at kt
91         !
92         !                                         ! ========================= !
93         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN      !    Add restoring term     !
94            !                                      ! ========================= !
95            !
96            IF( nn_sstr == 1 ) THEN                                   !* Temperature restoring term
97               DO_2D_11_11
98                  zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) * tmask(ji,jj,1)
99                  qns(ji,jj) = qns(ji,jj) + zqrp
100                  qrp(ji,jj) = zqrp
101               END_2D
102            ENDIF
103            !
104            IF( nn_sssr /= 0 .AND. nn_sssr_ice /= 1 ) THEN
105              ! use fraction of ice ( fr_i ) to adjust relaxation under ice if nn_sssr_ice .ne. 1
106              ! n.b. coefice is initialised and fixed to 1._wp if nn_sssr_ice = 1
107               DO_2D_11_11
108                  SELECT CASE ( nn_sssr_ice )
109                    CASE ( 0 )    ;  coefice(ji,jj) = 1._wp - fr_i(ji,jj)              ! no/reduced damping under ice
110                    CASE  DEFAULT ;  coefice(ji,jj) = 1._wp + ( nn_sssr_ice - 1 ) * fr_i(ji,jj) ! reinforced damping (x nn_sssr_ice) under ice )
111                  END SELECT
112               END_2D
113            ENDIF
114            !
115            IF( nn_sssr == 1 ) THEN                                   !* Salinity damping term (salt flux only (sfx))
116               zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
117               DO_2D_11_11
118                  zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
119                     &        *   coefice(ji,jj)            &      ! Optional control of damping under sea-ice
120                     &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) * tmask(ji,jj,1)
121                  sfx(ji,jj) = sfx(ji,jj) + zerp                 ! salt flux
122                  erp(ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20 ) ! converted into an equivalent volume flux (diagnostic only)
123               END_2D
124               !
125            ELSEIF( nn_sssr == 2 ) THEN                               !* Salinity damping term (volume flux (emp) and associated heat flux (qns)
126               zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
127               zerp_bnd = rn_sssr_bnd / rday                          !       -              -   
128               DO_2D_11_11
129                  zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
130                     &        *   coefice(ji,jj)            &      ! Optional control of damping under sea-ice
131                     &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )   &
132                     &        / MAX(  sss_m(ji,jj), 1.e-20   ) * tmask(ji,jj,1)
133                  IF( ln_sssr_bnd )   zerp = SIGN( 1., zerp ) * MIN( zerp_bnd, ABS(zerp) )
134                  emp(ji,jj) = emp (ji,jj) + zerp
135                  qns(ji,jj) = qns(ji,jj) - zerp * rcp * sst_m(ji,jj)
136                  erp(ji,jj) = zerp
137               END_2D
138            ENDIF
139            !
140         ENDIF
141         !
142      ENDIF
143      !
144   END SUBROUTINE sbc_ssr
145
146 
147   SUBROUTINE sbc_ssr_init
148      !!---------------------------------------------------------------------
149      !!                  ***  ROUTINE sbc_ssr_init  ***
150      !!
151      !! ** Purpose :   initialisation of surface damping term
152      !!
153      !! ** Method  : - Read namelist namsbc_ssr
154      !!              - Read observed SST and/or SSS if required
155      !!---------------------------------------------------------------------
156      INTEGER  ::   ji, jj   ! dummy loop indices
157      REAL(wp) ::   zerp     ! local scalar for evaporation damping
158      REAL(wp) ::   zqrp     ! local scalar for heat flux damping
159      REAL(wp) ::   zsrp     ! local scalar for unit conversion of rn_deds factor
160      REAL(wp) ::   zerp_bnd ! local scalar for unit conversion of rn_epr_max factor
161      INTEGER  ::   ierror   ! return error code
162      !!
163      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
164      TYPE(FLD_N) ::   sn_sst, sn_sss        ! informations about the fields to be read
165      NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, &
166              & sn_sss, ln_sssr_bnd, rn_sssr_bnd, nn_sssr_ice
167      INTEGER     ::  ios
168      !!----------------------------------------------------------------------
169      !
170      IF(lwp) THEN
171         WRITE(numout,*)
172         WRITE(numout,*) 'sbc_ssr : SST and/or SSS damping term '
173         WRITE(numout,*) '~~~~~~~ '
174      ENDIF
175      !
176      READ  ( numnam_ref, namsbc_ssr, IOSTAT = ios, ERR = 901)
177901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_ssr in reference namelist' )
178
179      READ  ( numnam_cfg, namsbc_ssr, IOSTAT = ios, ERR = 902 )
180902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_ssr in configuration namelist' )
181      IF(lwm) WRITE ( numond, namsbc_ssr )
182
183      IF(lwp) THEN                 !* control print
184         WRITE(numout,*) '   Namelist namsbc_ssr :'
185         WRITE(numout,*) '      SST restoring term (Yes=1)             nn_sstr        = ', nn_sstr
186         WRITE(numout,*) '         dQ/dT (restoring magnitude on SST)     rn_dqdt     = ', rn_dqdt, ' W/m2/K'
187         WRITE(numout,*) '      SSS damping term (Yes=1, salt   flux)  nn_sssr        = ', nn_sssr
188         WRITE(numout,*) '                       (Yes=2, volume flux) '
189         WRITE(numout,*) '         dE/dS (restoring magnitude on SST)     rn_deds     = ', rn_deds, ' mm/day'
190         WRITE(numout,*) '         flag to bound erp term                 ln_sssr_bnd = ', ln_sssr_bnd
191         WRITE(numout,*) '         ABS(Max./Min.) erp threshold           rn_sssr_bnd = ', rn_sssr_bnd, ' mm/day'
192         WRITE(numout,*) '      Cntrl of surface restoration under ice nn_sssr_ice    = ', nn_sssr_ice
193         WRITE(numout,*) '          ( 0 = no restoration under ice)'
194         WRITE(numout,*) '          ( 1 = restoration everywhere  )'
195         WRITE(numout,*) '          (>1 = enhanced restoration under ice  )'
196      ENDIF
197      !
198      IF( nn_sstr == 1 ) THEN      !* set sf_sst structure & allocate arrays
199         !
200         ALLOCATE( sf_sst(1), STAT=ierror )
201         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst structure' )
202         ALLOCATE( sf_sst(1)%fnow(jpi,jpj,1), STAT=ierror )
203         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst now array' )
204         !
205         ! fill sf_sst with sn_sst and control print
206         CALL fld_fill( sf_sst, (/ sn_sst /), cn_dir, 'sbc_ssr', 'SST restoring term toward SST data', 'namsbc_ssr', no_print )
207         IF( sf_sst(1)%ln_tint )   ALLOCATE( sf_sst(1)%fdta(jpi,jpj,1,2), STAT=ierror )
208         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst data array' )
209         !
210      ENDIF
211      !
212      IF( nn_sssr >= 1 ) THEN      !* set sf_sss structure & allocate arrays
213         !
214         ALLOCATE( sf_sss(1), STAT=ierror )
215         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss structure' )
216         ALLOCATE( sf_sss(1)%fnow(jpi,jpj,1), STAT=ierror )
217         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss now array' )
218         !
219         ! fill sf_sss with sn_sss and control print
220         CALL fld_fill( sf_sss, (/ sn_sss /), cn_dir, 'sbc_ssr', 'SSS restoring term toward SSS data', 'namsbc_ssr', no_print )
221         IF( sf_sss(1)%ln_tint )   ALLOCATE( sf_sss(1)%fdta(jpi,jpj,1,2), STAT=ierror )
222         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss data array' )
223         !
224      ENDIF
225      !
226      coefice(:,:) = 1._wp         !  Initialise coefice to 1._wp ; will not need to be changed if nn_sssr_ice=1
227      !                            !* Initialize qrp and erp if no restoring
228      IF( nn_sstr /= 1                   )   qrp(:,:) = 0._wp
229      IF( nn_sssr /= 1 .OR. nn_sssr /= 2 )   erp(:,:) = 0._wp
230      !
231   END SUBROUTINE sbc_ssr_init
232         
233   INTEGER FUNCTION sbc_ssr_alloc()
234      !!----------------------------------------------------------------------
235      !!               ***  FUNCTION sbc_ssr_alloc  ***
236      !!----------------------------------------------------------------------
237      sbc_ssr_alloc = 0       ! set to zero if no array to be allocated
238      IF( .NOT. ALLOCATED( erp ) ) THEN
239         ALLOCATE( qrp(jpi,jpj), erp(jpi,jpj), coefice(jpi,jpj), STAT= sbc_ssr_alloc )
240         !
241         IF( lk_mpp                  )   CALL mpp_sum ( 'sbcssr', sbc_ssr_alloc )
242         IF( sbc_ssr_alloc /= 0 )   CALL ctl_warn('sbc_ssr_alloc: failed to allocate arrays.')
243         !
244      ENDIF
245   END FUNCTION
246     
247   !!======================================================================
248END MODULE sbcssr
Note: See TracBrowser for help on using the repository browser.