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.
icesbc.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icesbc.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: 17.7 KB
Line 
1MODULE icesbc
2   !!======================================================================
3   !!                       ***  MODULE  icesbc  ***
4   !! Sea-Ice :   air-ice sbc fields
5   !!=====================================================================
6   !! History :  4.0  !  2017-08  (C. Rousset)       Original code
7   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
8   !!----------------------------------------------------------------------
9#if defined key_si3
10   !!----------------------------------------------------------------------
11   !!   'key_si3' :                                     SI3 sea-ice model
12   !!----------------------------------------------------------------------
13   USE oce            ! ocean dynamics and tracers
14   USE dom_oce        ! ocean space and time domain
15   USE ice            ! sea-ice: variables
16   USE sbc_oce        ! Surface boundary condition: ocean fields
17   USE sbc_ice        ! Surface boundary condition: ice   fields
18   USE usrdef_sbc     ! Surface boundary condition: user defined
19   USE sbcblk         ! Surface boundary condition: bulk
20   USE sbccpl         ! Surface boundary condition: coupled interface
21   USE icealb         ! sea-ice: albedo
22   !
23   USE in_out_manager ! I/O manager
24   USE iom            ! I/O manager library
25   USE lib_mpp        ! MPP library
26   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
27   USE lbclnk         ! lateral boundary conditions (or mpp links)
28   USE timing         ! Timing
29   USE fldread        !!GS: needed by agrif
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC ice_sbc_tau   ! called by icestp.F90
35   PUBLIC ice_sbc_flx   ! called by icestp.F90
36   PUBLIC ice_sbc_init  ! called by icestp.F90
37
38   !! * Substitutions
39#  include "vectopt_loop_substitute.h90"
40#  include "do_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
43   !! $Id$
44   !! Software governed by the CeCILL license (see ./LICENSE)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE ice_sbc_tau( kt, ksbc, utau_ice, vtau_ice )
49      !!-------------------------------------------------------------------
50      !!                  ***  ROUTINE ice_sbc_tau  ***
51      !!
52      !! ** Purpose : provide surface boundary condition for sea ice (momentum)
53      !!
54      !! ** Action  : It provides the following fields:
55      !!              utau_ice, vtau_ice : surface ice stress (U- & V-points) [N/m2]
56      !!-------------------------------------------------------------------
57      INTEGER                     , INTENT(in   ) ::   kt                   ! ocean time step
58      INTEGER                     , INTENT(in   ) ::   ksbc                 ! type of sbc flux
59      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   utau_ice, vtau_ice   ! air-ice stress   [N/m2]
60      !!
61      INTEGER  ::   ji, jj                 ! dummy loop index
62      REAL(wp), DIMENSION(jpi,jpj) ::   zutau_ice, zvtau_ice 
63      !!-------------------------------------------------------------------
64      !
65      IF( ln_timing )   CALL timing_start('ice_sbc')
66      !
67      IF( kt == nit000 .AND. lwp ) THEN
68         WRITE(numout,*)
69         WRITE(numout,*)'ice_sbc_tau: Surface boundary condition for sea ice (momentum)'
70         WRITE(numout,*)'~~~~~~~~~~~~~~~'
71      ENDIF
72      !
73      SELECT CASE( ksbc )
74         CASE( jp_usr     )   ;    CALL usrdef_sbc_ice_tau( kt )                 ! user defined formulation
75         CASE( jp_blk     )   ;    CALL blk_ice_1( sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),   &
76            &                                      sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1),   &
77            &                                      sf(jp_slp )%fnow(:,:,1), u_ice, v_ice, tm_su    ,   &   ! inputs
78            &                                      putaui = utau_ice, pvtaui = vtau_ice            )       ! outputs                             
79 !        CASE( jp_abl     )    utau_ice & vtau_ice are computed in ablmod
80         CASE( jp_purecpl )   ;    CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled      formulation
81      END SELECT
82      !
83      IF( ln_mixcpl) THEN                                                        ! Case of a mixed Bulk/Coupled formulation
84                                   CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
85         DO_2D_00_00
86            utau_ice(ji,jj) = utau_ice(ji,jj) * xcplmask(ji,jj,0) + zutau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
87            vtau_ice(ji,jj) = vtau_ice(ji,jj) * xcplmask(ji,jj,0) + zvtau_ice(ji,jj) * ( 1. - xcplmask(ji,jj,0) )
88         END_2D
89         CALL lbc_lnk_multi( 'icesbc', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
90      ENDIF
91      !
92      IF( ln_timing )   CALL timing_stop('ice_sbc')
93      !
94   END SUBROUTINE ice_sbc_tau
95
96   
97   SUBROUTINE ice_sbc_flx( kt, ksbc )
98      !!-------------------------------------------------------------------
99      !!                  ***  ROUTINE ice_sbc_flx  ***
100      !!
101      !! ** Purpose : provide surface boundary condition for sea ice (flux)
102      !!
103      !! ** Action  : It provides the following fields used in sea ice model:
104      !!                emp_oce , emp_ice                        = E-P over ocean and sea ice                    [Kg/m2/s]
105      !!                sprecip                                  = solid precipitation                           [Kg/m2/s]
106      !!                evap_ice                                 = sublimation                                   [Kg/m2/s]
107      !!                qsr_tot , qns_tot                        = solar & non solar heat flux (total)           [W/m2]
108      !!                qsr_ice , qns_ice                        = solar & non solar heat flux over ice          [W/m2]
109      !!                dqns_ice                                 = non solar  heat sensistivity                  [W/m2]
110      !!                qemp_oce, qemp_ice, qprec_ice, qevap_ice = sensible heat (associated with evap & precip) [W/m2]
111      !!            + some fields that are not used outside this module:
112      !!                qla_ice                                  = latent heat flux over ice                     [W/m2]
113      !!                dqla_ice                                 = latent heat sensistivity                      [W/m2]
114      !!                tprecip                                  = total  precipitation                          [Kg/m2/s]
115      !!                alb_ice                                  = albedo above sea ice
116      !!-------------------------------------------------------------------
117      INTEGER, INTENT(in) ::   kt     ! ocean time step
118      INTEGER, INTENT(in) ::   ksbc   ! flux formulation (user defined, bulk or Pure Coupled)
119      !
120      INTEGER  ::   ji, jj, jl      ! dummy loop index
121      REAL(wp) ::   zmiss_val       ! missing value retrieved from xios
122      REAL(wp), DIMENSION(jpi,jpj,jpl)              ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
123      REAL(wp), DIMENSION(:,:)        , ALLOCATABLE ::   zalb, zmsk00      ! 2D workspace
124      !!--------------------------------------------------------------------
125      !
126      IF( ln_timing )   CALL timing_start('ice_sbc_flx')
127
128      IF( kt == nit000 .AND. lwp ) THEN
129         WRITE(numout,*)
130         WRITE(numout,*)'ice_sbc_flx: Surface boundary condition for sea ice (flux)'
131         WRITE(numout,*)'~~~~~~~~~~~~~~~'
132      ENDIF
133
134      ! get missing value from xml
135      CALL iom_miss_val( "icetemp", zmiss_val )
136
137      ! --- cloud-sky and overcast-sky ice albedos --- !
138      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os )
139
140      ! albedo depends on cloud fraction because of non-linear spectral effects
141!!gm cldf_ice is a real, DOCTOR naming rule: start with cd means CHARACTER passed in argument !
142      alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
143      !
144      !
145      SELECT CASE( ksbc )   !== fluxes over sea ice ==!
146      !
147      CASE( jp_usr )              !--- user defined formulation
148                                  CALL usrdef_sbc_ice_flx( kt, h_s, h_i )
149      CASE( jp_blk, jp_abl )  !--- bulk formulation & ABL formulation
150                                  CALL blk_ice_2    ( t_su, h_s, h_i, alb_ice, sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1),    &
151            &                                           sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) )    !
152         IF( ln_mixcpl        )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
153         IF( nn_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
154         !                        !    compute conduction flux and surface temperature (as in Jules surface module)
155         IF( ln_cndflx .AND. .NOT.ln_cndemulate ) &
156            &                     CALL blk_ice_qcn    ( ln_virtual_itd, t_su, t_bo, h_s, h_i )
157      CASE ( jp_purecpl )         !--- coupled formulation
158                                  CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
159         IF( nn_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
160      END SELECT
161
162      !--- output ice albedo and surface albedo ---!
163      IF( iom_use('icealb') .OR. iom_use('albedo') ) THEN
164
165         ALLOCATE( zalb(jpi,jpj), zmsk00(jpi,jpj) )
166
167         WHERE( at_i_b < 1.e-03 )
168            zmsk00(:,:) = 0._wp
169            zalb  (:,:) = rn_alb_oce
170         ELSEWHERE
171            zmsk00(:,:) = 1._wp           
172            zalb  (:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b
173         END WHERE
174         ! ice albedo
175         CALL iom_put( 'icealb' , zalb * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )
176         ! ice+ocean albedo
177         zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b )
178         CALL iom_put( 'albedo' , zalb )
179
180         DEALLOCATE( zalb, zmsk00 )
181
182      ENDIF
183      !
184      IF( ln_timing )   CALL timing_stop('ice_sbc_flx')
185      !
186   END SUBROUTINE ice_sbc_flx
187
188
189   SUBROUTINE ice_flx_dist( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_flxdist )
190      !!-------------------------------------------------------------------
191      !!                  ***  ROUTINE ice_flx_dist  ***
192      !!
193      !! ** Purpose :   update the ice surface boundary condition by averaging
194      !!              and/or redistributing fluxes on ice categories
195      !!
196      !! ** Method  :   average then redistribute
197      !!
198      !! ** Action  :   depends on k_flxdist
199      !!                = -1  Do nothing (needs N(cat) fluxes)
200      !!                =  0  Average N(cat) fluxes then apply the average over the N(cat) ice
201      !!                =  1  Average N(cat) fluxes then redistribute over the N(cat) ice
202      !!                                                 using T-ice and albedo sensitivity
203      !!                =  2  Redistribute a single flux over categories
204      !!-------------------------------------------------------------------
205      INTEGER                   , INTENT(in   ) ::   k_flxdist  ! redistributor
206      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptn_ice    ! ice surface temperature
207      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb_ice   ! ice albedo
208      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqns_ice   ! non solar flux
209      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqsr_ice   ! net solar flux
210      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdqn_ice   ! non solar flux sensitivity
211      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pevap_ice  ! sublimation
212      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdevap_ice ! sublimation sensitivity
213      !
214      INTEGER  ::   jl      ! dummy loop index
215      !
216      REAL(wp), DIMENSION(jpi,jpj) ::   z1_at_i   ! inverse of concentration
217      !
218      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_qsr_m   ! Mean solar heat flux over all categories
219      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_qns_m   ! Mean non solar heat flux over all categories
220      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_evap_m  ! Mean sublimation over all categories
221      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_dqn_m   ! Mean d(qns)/dT over all categories
222      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_devap_m ! Mean d(evap)/dT over all categories
223      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zalb_m    ! Mean albedo over all categories
224      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztem_m    ! Mean temperature over all categories
225      !!----------------------------------------------------------------------
226      !
227      WHERE ( at_i (:,:) > 0._wp )   ; z1_at_i(:,:) = 1._wp / at_i (:,:)
228      ELSEWHERE                      ; z1_at_i(:,:) = 0._wp
229      END WHERE
230     
231      SELECT CASE( k_flxdist )       !==  averaged on all ice categories  ==!
232      !
233      CASE( 0 , 1 )
234         !
235         ALLOCATE( z_qns_m(jpi,jpj), z_qsr_m(jpi,jpj), z_dqn_m(jpi,jpj), z_evap_m(jpi,jpj), z_devap_m(jpi,jpj) ) 
236         !
237         z_qns_m  (:,:) = SUM( a_i(:,:,:) * pqns_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
238         z_qsr_m  (:,:) = SUM( a_i(:,:,:) * pqsr_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
239         z_dqn_m  (:,:) = SUM( a_i(:,:,:) * pdqn_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
240         z_evap_m (:,:) = SUM( a_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
241         z_devap_m(:,:) = SUM( a_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
242         DO jl = 1, jpl
243            pqns_ice  (:,:,jl) = z_qns_m (:,:)
244            pqsr_ice  (:,:,jl) = z_qsr_m (:,:)
245            pdqn_ice  (:,:,jl) = z_dqn_m  (:,:)
246            pevap_ice (:,:,jl) = z_evap_m(:,:)
247            pdevap_ice(:,:,jl) = z_devap_m(:,:)
248         END DO
249         !
250         DEALLOCATE( z_qns_m, z_qsr_m, z_dqn_m, z_evap_m, z_devap_m ) 
251         !
252      END SELECT
253      !
254      SELECT CASE( k_flxdist )       !==  redistribution on all ice categories  ==!
255      !
256      CASE( 1 , 2 )
257         !
258         ALLOCATE( zalb_m(jpi,jpj), ztem_m(jpi,jpj) ) 
259         !
260         zalb_m(:,:) = SUM( a_i(:,:,:) * palb_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
261         ztem_m(:,:) = SUM( a_i(:,:,:) * ptn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
262         DO jl = 1, jpl
263            pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice  (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
264            pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
265            pqsr_ice (:,:,jl) = pqsr_ice (:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )
266         END DO
267         !
268         DEALLOCATE( zalb_m, ztem_m ) 
269         !
270      END SELECT
271      !
272   END SUBROUTINE ice_flx_dist
273
274
275   SUBROUTINE ice_sbc_init
276      !!-------------------------------------------------------------------
277      !!                  ***  ROUTINE ice_sbc_init  ***
278      !!
279      !! ** Purpose :   Physical constants and parameters linked to the ice dynamics
280      !!     
281      !! ** Method  :   Read the namsbc namelist and check the ice-dynamic
282      !!              parameter values called at the first timestep (nit000)
283      !!
284      !! ** input   :   Namelist namsbc
285      !!-------------------------------------------------------------------
286      INTEGER ::   ios, ioptio   ! Local integer
287      !!
288      NAMELIST/namsbc/ rn_cio, rn_blow_s, nn_flxdist, ln_cndflx, ln_cndemulate
289      !!-------------------------------------------------------------------
290      !
291      READ  ( numnam_ice_ref, namsbc, IOSTAT = ios, ERR = 901)
292901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc in reference namelist' )
293      READ  ( numnam_ice_cfg, namsbc, IOSTAT = ios, ERR = 902 )
294902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc in configuration namelist' )
295      IF(lwm) WRITE( numoni, namsbc )
296      !
297      IF(lwp) THEN                     ! control print
298         WRITE(numout,*)
299         WRITE(numout,*) 'ice_sbc_init: ice parameters for ice dynamics '
300         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
301         WRITE(numout,*) '   Namelist namsbc:'
302         WRITE(numout,*) '      drag coefficient for oceanic stress              rn_cio        = ', rn_cio
303         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_blow_s     = ', rn_blow_s
304         WRITE(numout,*) '      Multicategory heat flux formulation              nn_flxdist    = ', nn_flxdist
305         WRITE(numout,*) '      Use conduction flux as surface condition         ln_cndflx     = ', ln_cndflx
306         WRITE(numout,*) '         emulate conduction flux                       ln_cndemulate = ', ln_cndemulate
307      ENDIF
308      !
309      IF(lwp) WRITE(numout,*)
310      SELECT CASE( nn_flxdist )         ! SI3 Multi-category heat flux formulation
311      CASE( -1  )
312         IF(lwp) WRITE(numout,*) '   SI3: use per-category fluxes (nn_flxdist = -1) '
313      CASE(  0  )
314         IF(lwp) WRITE(numout,*) '   SI3: use average per-category fluxes (nn_flxdist = 0) '
315      CASE(  1  )
316         IF(lwp) WRITE(numout,*) '   SI3: use average then redistribute per-category fluxes (nn_flxdist = 1) '
317         IF( ln_cpl )         CALL ctl_stop( 'ice_thd_init: the chosen nn_flxdist for SI3 in coupled mode must be /=1' )
318      CASE(  2  )
319         IF(lwp) WRITE(numout,*) '   SI3: Redistribute a single flux over categories (nn_flxdist = 2) '
320         IF( .NOT. ln_cpl )   CALL ctl_stop( 'ice_thd_init: the chosen nn_flxdist for SI3 in forced mode must be /=2' )
321      CASE DEFAULT
322         CALL ctl_stop( 'ice_thd_init: SI3 option, nn_flxdist, should be between -1 and 2' )
323      END SELECT
324      !
325   END SUBROUTINE ice_sbc_init
326
327#else
328   !!----------------------------------------------------------------------
329   !!   Default option :         Empty module         NO SI3 sea-ice model
330   !!----------------------------------------------------------------------
331#endif
332
333   !!======================================================================
334END MODULE icesbc
Note: See TracBrowser for help on using the repository browser.