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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF/isfcavmlt.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.

File size: 16.8 KB
Line 
1MODULE isfcavmlt
2   !!======================================================================
3   !!                       ***  MODULE  isfcavmlt  ***
4   !! ice shelf module :  update surface ocean boundary condition under ice
5   !!                   shelves
6   !!======================================================================
7   !! History :  4.0  !  2019-09  (P. Mathiot) Original code
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   isfcav_mlt    : compute or read ice shelf fwf/heat fluxes in the ice shelf cavity
12   !!----------------------------------------------------------------------
13
14   USE isf_oce                  ! ice shelf
15   USE isftbl , ONLY: isf_tbl   ! ice shelf depth average
16   USE isfutils,ONLY: debug     ! debug subroutine
17
18   USE dom_oce                            ! ocean space and time domain
19   USE phycst , ONLY: rcp, rau0, rau0_rcp ! physical constants
20   USE eosbn2 , ONLY: eos_fzp             ! equation of state
21
22   USE in_out_manager              ! I/O manager
23   USE iom        , ONLY: iom_put  ! I/O library
24   USE fldread    , ONLY: fld_read, FLD, FLD_N !
25   USE lib_fortran, ONLY: glob_sum !
26   USE lib_mpp    , ONLY: ctl_stop !
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   isfcav_mlt
32
33   !! * Substitutions
34#  include "do_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
37   !! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
38   !! Software governed by the CeCILL license (see ./LICENSE)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42! -------------------------------------------------------------------------------------------------------
43! -------------------------------- PUBLIC SUBROUTINE ----------------------------------------------------
44! -------------------------------------------------------------------------------------------------------
45
46   SUBROUTINE isfcav_mlt(kt, pgt, pgs , pttbl, pstbl, &
47      &                           pqhc, pqoce, pqfwf  )
48      !!----------------------------------------------------------------------
49      !!
50      !!                          ***  ROUTINE isfcav_mlt  ***
51      !!
52      !! ** Purpose    : compute or read ice shelf fwf/heat fluxes in the ice shelf cavity
53      !!
54      !!-------------------------- OUT -------------------------------------
55      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! heat and fwf fluxes
56      !!-------------------------- IN  -------------------------------------
57      INTEGER, INTENT(in) :: kt
58      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt  , pgs    ! gamma t and gamma s
59      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl  ! top boundary layer tracer
60      !!---------------------------------------------------------------------
61      !
62      ! compute latent heat and melt (2d)
63      SELECT CASE ( cn_isfcav_mlt )
64      CASE ( 'spe' )   ! ice shelf melt specified (read input file, and heat fluxes derived from
65         CALL isfcav_mlt_spe( kt, pstbl,               &
66            &                  pqhc, pqoce, pqfwf  )
67      CASE ( '2eq' )   !  ISOMIP  formulation (2 equations) for volume flux (Hunter et al., 2006)
68         CALL isfcav_mlt_2eq( pgt, pttbl, pstbl,       &
69            &                  pqhc , pqoce, pqfwf )
70      CASE ( '3eq' )   ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
71         CALL isfcav_mlt_3eq( pgt, pgs , pttbl, pstbl, &
72            &                  pqhc, pqoce, pqfwf  )
73      CASE ( 'oasis' ) ! fwf pass trough oasis
74         CALL isfcav_mlt_oasis( kt, pstbl,              &
75            &                   pqhc, pqoce, pqfwf  )
76      CASE DEFAULT
77         CALL ctl_stop('STOP', 'unknown isf melt formulation : cn_isfcav (should not see this)')
78      END SELECT
79      !
80      IF (ln_isfdebug) THEN
81         CALL debug( 'isfcav_mlt qhc  :', pqhc (:,:) )
82         CALL debug( 'isfcav_mlt qoce :', pqoce(:,:) )
83         CALL debug( 'isfcav_mlt qfwf :', pqfwf(:,:) )
84      END IF
85      !
86   END SUBROUTINE isfcav_mlt
87
88! -------------------------------------------------------------------------------------------------------
89! -------------------------------- PRIVATE SUBROUTINE ---------------------------------------------------
90! -------------------------------------------------------------------------------------------------------
91
92   SUBROUTINE isfcav_mlt_spe(kt, pstbl,          &  ! <<== in
93      &                      pqhc , pqoce, pqfwf )  ! ==>> out
94      !!----------------------------------------------------------------------
95      !!
96      !!                          ***  ROUTINE isfcav_mlt_spe  ***
97      !!
98      !! ** Purpose    : - read ice shelf melt from forcing file
99      !!                 - compute ocea-ice heat flux (assuming it is equal to latent heat)
100      !!                 - compute heat content flux
101      !!---------------------------------------------------------------------
102      !!-------------------------- OUT -------------------------------------
103      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! heat content, latent heat and fwf fluxes
104      !!-------------------------- IN  -------------------------------------
105      INTEGER                     , INTENT(in   ) :: kt                  ! current time step
106      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pstbl               ! salinity in tbl
107      !!--------------------------------------------------------------------
108      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz                              ! tbl freezing temperature
109      !!--------------------------------------------------------------------
110      !
111      ! Compute freezing temperature
112      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
113      !
114      ! read input file
115      CALL fld_read ( kt, 1, sf_isfcav_fwf )
116      !
117      ! define fwf and qoce
118      ! ocean heat flux is assume to be equal to the latent heat
119      pqfwf(:,:) = - sf_isfcav_fwf(1)%fnow(:,:,1)      ! fwf                ( >0 out)
120      pqoce(:,:) = - pqfwf(:,:) * rLfusisf             ! ocean heat flux    ( >0 out)
121      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux  ( >0 out)
122      !
123      ! output freezing point at the interface
124      CALL iom_put('isftfrz_cav', ztfrz )
125      !
126   END SUBROUTINE isfcav_mlt_spe
127
128   SUBROUTINE isfcav_mlt_2eq(pgt , pttbl, pstbl, &  ! <<== in
129      &                      pqhc, pqoce, pqfwf  )  ! ==>> out
130      !!----------------------------------------------------------------------
131      !!
132      !!                          ***  ROUTINE isfcav_mlt_2eq  ***
133      !!
134      !! ** Purpose    : Compute ice shelf fwf/heqt fluxes using ISOMIP formulation (Hunter et al., 2006)
135      !!
136      !! ** Method     : The ice shelf melt latent heat is defined as being equal to the ocean/ice heat flux.
137      !!                 From this we can derived the fwf, ocean/ice heat flux and the heat content flux as being :
138      !!                   qfwf  = Gammat * Rau0 * Cp * ( Tw - Tfrz ) / Lf
139      !!                   qhoce = qlat
140      !!                   qhc   = qfwf * Cp * Tfrz
141      !!
142      !! ** Reference  : Hunter,  J.  R.:  Specification  for  test  models  of  ice  shelf  cavities, 
143      !!                 Tech.  Rep.  June,  Antarctic  Climate  &  Ecosystems  Cooperative  Research  Centre,  available  at: 
144      !!                 http://staff.acecrc.org.au/~bkgalton/ISOMIP/test_cavities.pdf (last access: 21 July 2016), 2006.
145      !!
146      !!-------------------------- OUT -------------------------------------
147      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! hean content, ocean-ice heat and fwf fluxes
148      !!-------------------------- IN  -------------------------------------
149      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt           ! temperature exchange coeficient
150      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl  ! temperature and salinity in top boundary layer
151      !!--------------------------------------------------------------------
152      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz         ! freezing temperature
153      REAL(wp), DIMENSION(jpi,jpj) :: zthd          ! thermal driving
154      !!--------------------------------------------------------------------
155      !
156      ! Calculate freezing temperature
157      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
158      !
159      ! thermal driving
160      zthd (:,:) = ( pttbl(:,:) - ztfrz(:,:) ) * mskisf_cav(:,:)
161      !
162      ! compute ocean-ice heat flux and then derive fwf assuming that ocean heat flux equal latent heat
163      pqfwf(:,:) = - pgt(:,:) * rau0_rcp * zthd(:,:) / rLfusisf    ! fresh water flux  ( > 0 out )
164      pqoce(:,:) = - pqfwf(:,:) * rLfusisf                         ! ocea-ice flux     ( > 0 out )
165      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp                 ! heat content flux ( > 0 out )
166      !
167      ! output thermal driving and freezinpoint at the ice shelf interface
168      CALL iom_put('isfthermald_cav', zthd )
169      CALL iom_put('isftfrz_cav'    , ztfrz )
170      !
171   END SUBROUTINE isfcav_mlt_2eq
172
173   SUBROUTINE isfcav_mlt_3eq(pgt, pgs , pttbl, pstbl, &  ! <<== in
174      &                           pqhc, pqoce, pqfwf  )  ! ==>> out
175      !!----------------------------------------------------------------------
176      !!
177      !!                          ***  ROUTINE isfcav_mlt_3eq  ***
178      !!
179      !! ** Purpose    : Compute ice shelf fwf/heqt fluxes using the 3 equation formulation
180      !!
181      !! ** Method     : The melt rate is determined considering the heat balance, the salt balance
182      !!                 at the phase change interface and a linearisation of the equation of state.
183      !!
184      !! ** Reference  : - Holland, D. M. and Jenkins, A.,
185      !!                   Modeling Thermodynamic Ice-Ocean Interactions at the Base of an Ice Shelf,
186      !!                   J. Phys. Oceanogr., 29, 1999.
187      !!                 - Asay-Davis, X. S., Cornford, S. L., Durand, G., Galton-Fenzi, B. K., Gladstone,
188      !!                   R. M., Gudmundsson, G. H., Hattermann, T., Holland, D. M., Holland, D., Holland,
189      !!                   P. R., Martin, D. F., Mathiot, P., Pattyn, F., and Seroussi, H.:
190      !!                   Experimental design for three interrelated marine ice sheet and ocean model intercomparison projects:
191      !!                   MISMIP v. 3 (MISMIP +), ISOMIP v. 2 (ISOMIP +) and MISOMIP v. 1 (MISOMIP1),
192      !!                   Geosci. Model Dev., 9, 2471-2497, https://doi.org/10.5194/gmd-9-2471-2016, 2016.
193      !!
194      !!-------------------------- OUT -------------------------------------
195      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! latent heat and fwf fluxes
196      !!-------------------------- IN  -------------------------------------
197      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt  , pgs          ! heat/salt exchange coeficient
198      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl        ! mean temperature and salinity in top boundary layer
199      !!--------------------------------------------------------------------
200      REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7       ! dummy local scalar for quadratic equation resolution
201      REAL(wp) :: zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac     ! dummy local scalar for quadratic equation resolution
202      REAL(wp) :: zeps = 1.e-20
203      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz         ! freezing point
204      REAL(wp), DIMENSION(jpi,jpj) :: zqcon         ! conductive flux through the ice shelf
205      REAL(wp), DIMENSION(jpi,jpj) :: zthd          ! thermal driving
206      !
207      INTEGER  ::   ji, jj     ! dummy loop indices
208      !!--------------------------------------------------------------------
209      !
210      ! compute upward heat flux zhtflx and upward water flux zwflx
211      ! Resolution of a 3d equation from equation 24, 25 and 26 (note conduction through the ice has been added to Eq 24)
212      DO_2D_11_11
213         !
214         ! compute coeficient to solve the 2nd order equation
215         zeps1 = rau0_rcp * pgt(ji,jj)
216         zeps2 = rLfusisf * rau0 * pgs(ji,jj)
217         zeps3 = rhoisf * rcpisf * rkappa / MAX(risfdep(ji,jj),zeps)
218         zeps4 = risf_lamb2 + risf_lamb3 * risfdep(ji,jj)
219         zeps6 = zeps4 - pttbl(ji,jj)
220         zeps7 = zeps4 - rtsurf
221         !
222         ! solve the 2nd order equation to find zsfrz
223         zaqe  = risf_lamb1 * (zeps1 + zeps3)
224         zaqer = 0.5_wp / MIN(zaqe,-zeps)
225         zbqe  = zeps1 * zeps6 + zeps3 * zeps7 - zeps2
226         zcqe  = zeps2 * pstbl(ji,jj)
227         zdis  = zbqe * zbqe - 4.0_wp * zaqe * zcqe               
228         !
229         ! Presumably zdis can never be negative because gammas is very small compared to gammat
230         zsfrz=(-zbqe - SQRT(zdis)) * zaqer
231         IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe + SQRT(zdis)) * zaqer  ! check this if this if is needed
232         !
233         ! compute t freeze (eq. 25)
234         ztfrz(ji,jj) = zeps4 + risf_lamb1 * zsfrz
235         !
236         ! thermal driving
237         zthd(ji,jj) = ( pttbl(ji,jj) - ztfrz(ji,jj) )
238         !
239         ! compute the upward water and heat flux (eq. 24 and eq. 26)
240         pqfwf(ji,jj) = rau0     * pgs(ji,jj) * ( zsfrz - pstbl(ji,jj) ) / MAX(zsfrz,zeps) ! fresh water flux    (> 0 out)
241         pqoce(ji,jj) = rau0_rcp * pgt(ji,jj) * zthd (ji,jj)                               ! ocean-ice heat flux (> 0 out)
242         pqhc (ji,jj) = rcp      * pqfwf(ji,jj) * ztfrz(ji,jj)                             ! heat content   flux (> 0 out)
243         !
244         zqcon(ji,jj) = zeps3 * ( ztfrz(ji,jj) - rtsurf )
245         !
246      END_2D
247      !
248      ! output conductive heat flux through the ice
249      CALL iom_put('qconisf', zqcon(:,:) * mskisf_cav(:,:) )
250      !
251      ! output thermal driving and freezing point at the interface
252      CALL iom_put('isfthermald_cav', zthd (:,:) * mskisf_cav(:,:) )
253      CALL iom_put('isftfrz_cav'    , ztfrz(:,:) * mskisf_cav(:,:) )
254      !
255   END SUBROUTINE isfcav_mlt_3eq
256
257   SUBROUTINE isfcav_mlt_oasis(kt, pstbl,          &  ! <<== in
258      &                        pqhc , pqoce, pqfwf )  ! ==>> out
259      !!----------------------------------------------------------------------
260      !!                          ***  ROUTINE isfcav_mlt_oasis  ***
261      !!
262      !! ** Purpose    : scale the fwf read from input file by the total amount received by the sbccpl interface
263      !!
264      !! ** Purpose    : - read ice shelf melt from forcing file => pattern
265      !!                 - total amount of fwf is given by sbccpl (fwfisf_cpl)
266      !!                 - scale fwf and compute heat fluxes
267      !!
268      !!---------------------------------------------------------------------
269      !!-------------------------- OUT -------------------------------------
270      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! heat content, latent heat and fwf fluxes
271      !!-------------------------- IN  -------------------------------------
272      INTEGER                     , INTENT(in   ) :: kt                  ! current time step
273      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pstbl               ! salinity in tbl
274      !!--------------------------------------------------------------------
275      REAL(wp)                     :: zfwf_fld, zfwf_oasis               ! total fwf in the forcing fields (pattern) and from the oasis interface (amount)
276      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz                              ! tbl freezing temperature
277      REAL(wp), DIMENSION(jpi,jpj) :: zfwf                               ! 2d fwf map after scaling
278      !!--------------------------------------------------------------------
279      !
280      ! Calculate freezing temperature
281      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
282      !
283      ! read input file
284      CALL fld_read ( kt, 1, sf_isfcav_fwf )
285      !
286      ! ice shelf 2d map
287      zfwf(:,:) = - sf_isfcav_fwf(1)%fnow(:,:,1)
288      !
289      ! compute glob sum from input file
290      ! (PM) should consider delay sum as in fwb (1 time step offset if I well understood)
291      zfwf_fld = glob_sum('isfcav_mlt', e1e2t(:,:) * zfwf(:,:))
292      !
293      ! compute glob sum from atm->oce ice shelf fwf
294      ! (PM) should consider delay sum as in fwb (1 time step offset if I well understood)
295      zfwf_oasis = glob_sum('isfcav_mlt', e1e2t(:,:) * fwfisf_oasis(:,:))
296      !
297      ! scale fwf
298      zfwf(:,:) = zfwf(:,:) * zfwf_oasis / zfwf_fld
299      !
300      ! define fwf and qoce
301      ! ocean heat flux is assume to be equal to the latent heat
302      pqfwf(:,:) =   zfwf(:,:)                         ! fwf                ( >0 out)
303      pqoce(:,:) = - pqfwf(:,:) * rLfusisf             ! ocean heat flux    ( >0 out)
304      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux  ( >0 out)
305      !
306      CALL iom_put('isftfrz_cav', ztfrz * mskisf_cav(:,:) )
307      !
308   END SUBROUTINE isfcav_mlt_oasis
309
310END MODULE isfcavmlt
Note: See TracBrowser for help on using the repository browser.