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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/iceupdate.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: 26.6 KB
Line 
1MODULE iceupdate
2   !!======================================================================
3   !!                       ***  MODULE iceupdate   ***
4   !!  Sea-ice :   computation of the flux at the sea ice/ocean interface
5   !!======================================================================
6   !! History :  4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
7   !!----------------------------------------------------------------------
8#if defined key_si3
9   !!----------------------------------------------------------------------
10   !!   'key_si3'                                       SI3 sea-ice model
11   !!----------------------------------------------------------------------
12   !!   ice_update_alloc : allocate the iceupdate arrays
13   !!   ice_update_init  : initialisation
14   !!   ice_update_flx   : updates mass, heat and salt fluxes at the ocean surface
15   !!   ice_update_tau   : update i- and j-stresses, and its modulus at the ocean surface
16   !!----------------------------------------------------------------------
17   USE phycst         ! physical constants
18   USE dom_oce        ! ocean domain
19   USE ice            ! sea-ice: variables
20   USE sbc_ice        ! Surface boundary condition: ice   fields
21   USE sbc_oce        ! Surface boundary condition: ocean fields
22   USE sbccpl         ! Surface boundary condition: coupled interface
23   USE icealb         ! sea-ice: albedo parameters
24   USE traqsr         ! add penetration of solar flux in the calculation of heat budget
25   USE icectl         ! sea-ice: control prints
26   USE bdy_oce , ONLY : ln_bdy
27   !
28   USE in_out_manager ! I/O manager
29   USE iom            ! I/O manager library
30   USE lib_mpp        ! MPP library
31   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
32   USE lbclnk         ! lateral boundary conditions (or mpp links)
33   USE timing         ! Timing
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   ice_update_init   ! called by ice_init
39   PUBLIC   ice_update_flx    ! called by ice_stp
40   PUBLIC   ice_update_tau    ! called by ice_stp
41
42   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress     [N/m2]
43   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmod_io              ! modulus of the ice-ocean velocity   [m/s]
44
45   !! * Substitutions
46#  include "vectopt_loop_substitute.h90"
47#  include "do_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
50   !! $Id$
51   !! Software governed by the CeCILL license (see ./LICENSE)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   INTEGER FUNCTION ice_update_alloc()
56      !!-------------------------------------------------------------------
57      !!             ***  ROUTINE ice_update_alloc ***
58      !!-------------------------------------------------------------------
59      ALLOCATE( utau_oce(jpi,jpj), vtau_oce(jpi,jpj), tmod_io(jpi,jpj), STAT=ice_update_alloc )
60      !
61      CALL mpp_sum( 'iceupdate', ice_update_alloc )
62      IF( ice_update_alloc /= 0 )   CALL ctl_stop( 'STOP', 'ice_update_alloc: failed to allocate arrays' )
63      !
64   END FUNCTION ice_update_alloc
65
66
67   SUBROUTINE ice_update_flx( kt )
68      !!-------------------------------------------------------------------
69      !!                ***  ROUTINE ice_update_flx ***
70      !! 
71      !! ** Purpose :   Update the surface ocean boundary condition for heat
72      !!                salt and mass over areas where sea-ice is non-zero
73      !!         
74      !! ** Action  : - computes the heat and freshwater/salt fluxes
75      !!                at the ice-ocean interface.
76      !!              - Update the ocean sbc
77      !!     
78      !! ** Outputs : - qsr     : sea heat flux:     solar
79      !!              - qns     : sea heat flux: non solar
80      !!              - emp     : freshwater budget: volume flux
81      !!              - sfx     : salt flux
82      !!              - fr_i    : ice fraction
83      !!              - tn_ice  : sea-ice surface temperature
84      !!              - alb_ice : sea-ice albedo (recomputed only for coupled mode)
85      !!
86      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
87      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
88      !!              These refs are now obsolete since everything has been revised
89      !!              The ref should be Rousset et al., 2015
90      !!---------------------------------------------------------------------
91      INTEGER, INTENT(in) ::   kt   ! number of iteration
92      !
93      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
94      REAL(wp) ::   zqmass           ! Heat flux associated with mass exchange ice->ocean (W.m-2)
95      REAL(wp) ::   zqsr             ! New solar flux received by the ocean
96      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d                  ! 2D workspace
97      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_cs, zalb_os     ! 3D workspace
98      !!---------------------------------------------------------------------
99      IF( ln_timing )   CALL timing_start('ice_update')
100
101      IF( kt == nit000 .AND. lwp ) THEN
102         WRITE(numout,*)
103         WRITE(numout,*)'ice_update_flx: update fluxes (mass, salt and heat) at the ice-ocean interface'
104         WRITE(numout,*)'~~~~~~~~~~~~~~'
105      ENDIF
106
107      ! --- case we bypass ice thermodynamics --- !
108      IF( .NOT. ln_icethd ) THEN   ! we suppose ice is impermeable => ocean is isolated from atmosphere
109         qt_atm_oi  (:,:)   = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
110         qt_oce_ai  (:,:)   = ( 1._wp - at_i_b(:,:) ) *   qns_oce(:,:)                  + qemp_oce(:,:)
111         emp_ice    (:,:)   = 0._wp
112         qemp_ice   (:,:)   = 0._wp
113         qevap_ice  (:,:,:) = 0._wp
114      ENDIF
115     
116      DO_2D_11_11
117
118         ! Solar heat flux reaching the ocean = zqsr (W.m-2)
119         !---------------------------------------------------
120         zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - qtr_ice_bot(ji,jj,:) ) )
121
122         ! Total heat flux reaching the ocean = qt_oce_ai (W.m-2)
123         !---------------------------------------------------
124         zqmass           = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC)
125         qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + zqmass + zqsr
126
127         ! Add the residual from heat diffusion equation and sublimation (W.m-2)
128         !----------------------------------------------------------------------
129         qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + hfx_err_dif(ji,jj) +   &
130            &             ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) )
131
132         ! New qsr and qns used to compute the oceanic heat flux at the next time step
133         !----------------------------------------------------------------------------
134         qsr(ji,jj) = zqsr                                     
135         qns(ji,jj) = qt_oce_ai(ji,jj) - zqsr             
136
137         ! Mass flux at the atm. surface       
138         !-----------------------------------
139         wfx_sub(ji,jj) = wfx_snw_sub(ji,jj) + wfx_ice_sub(ji,jj)
140
141         ! Mass flux at the ocean surface     
142         !------------------------------------
143         !  case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED)
144         !  -------------------------------------------------------------------------------------
145         !  The idea of this approach is that the system that we consider is the ICE-OCEAN system
146         !  Thus  FW  flux  =  External ( E-P+snow melt)
147         !       Salt flux  =  Exchanges in the ice-ocean system then converted into FW
148         !                     Associated to Ice formation AND Ice melting
149         !                     Even if i see Ice melting as a FW and SALT flux
150         !       
151         ! mass flux from ice/ocean
152         wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   &
153            &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj)
154
155         ! add the snow melt water to snow mass flux to the ocean
156         wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj)
157
158         ! mass flux at the ocean/ice interface
159         fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) )              ! F/M mass flux save at least for biogeochemical model
160         emp(ji,jj)    = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange)
161
162
163         ! Salt flux at the ocean surface     
164         !------------------------------------------
165         sfx(ji,jj) = sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj)   &
166            &       + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj)
167         
168         ! Mass of snow and ice per unit area   
169         !----------------------------------------
170         snwice_mass_b(ji,jj) = snwice_mass(ji,jj)       ! save mass from the previous ice time step
171         !                                               ! new mass per unit area
172         snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj)  ) 
173         !                                               ! time evolution of snow+ice mass
174         snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_rdtice
175         
176      END_2D
177
178      ! Storing the transmitted variables
179      !----------------------------------
180      fr_i  (:,:)   = at_i(:,:)             ! Sea-ice fraction           
181      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                     
182
183      ! Snow/ice albedo (only if sent to coupler, useless in forced mode)
184      !------------------------------------------------------------------
185      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
186      !
187      alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
188      !
189      IF( lrst_ice ) THEN                       !* write snwice_mass fields in the restart file
190         CALL update_rst( 'WRITE', kt )
191      ENDIF
192      !
193      ! output all fluxes
194      !------------------
195      !
196      ! --- salt fluxes [kg/m2/s] --- !
197      !                           ! sfxice =  sfxbog + sfxbom + sfxsum + sfxsni + sfxopw + sfxres + sfxdyn + sfxbri + sfxsub + sfxlam
198      IF( iom_use('sfxice'  ) )   CALL iom_put( 'sfxice', sfx     * 1.e-03 )   ! salt flux from total ice growth/melt
199      IF( iom_use('sfxbog'  ) )   CALL iom_put( 'sfxbog', sfx_bog * 1.e-03 )   ! salt flux from bottom growth
200      IF( iom_use('sfxbom'  ) )   CALL iom_put( 'sfxbom', sfx_bom * 1.e-03 )   ! salt flux from bottom melting
201      IF( iom_use('sfxsum'  ) )   CALL iom_put( 'sfxsum', sfx_sum * 1.e-03 )   ! salt flux from surface melting
202      IF( iom_use('sfxlam'  ) )   CALL iom_put( 'sfxlam', sfx_lam * 1.e-03 )   ! salt flux from lateral melting
203      IF( iom_use('sfxsni'  ) )   CALL iom_put( 'sfxsni', sfx_sni * 1.e-03 )   ! salt flux from snow ice formation
204      IF( iom_use('sfxopw'  ) )   CALL iom_put( 'sfxopw', sfx_opw * 1.e-03 )   ! salt flux from open water formation
205      IF( iom_use('sfxdyn'  ) )   CALL iom_put( 'sfxdyn', sfx_dyn * 1.e-03 )   ! salt flux from ridging rafting
206      IF( iom_use('sfxbri'  ) )   CALL iom_put( 'sfxbri', sfx_bri * 1.e-03 )   ! salt flux from brines
207      IF( iom_use('sfxres'  ) )   CALL iom_put( 'sfxres', sfx_res * 1.e-03 )   ! salt flux from undiagnosed processes
208      IF( iom_use('sfxsub'  ) )   CALL iom_put( 'sfxsub', sfx_sub * 1.e-03 )   ! salt flux from sublimation
209
210      ! --- mass fluxes [kg/m2/s] --- !
211      CALL iom_put( 'emp_oce', emp_oce )   ! emp over ocean (taking into account the snow blown away from the ice)
212      CALL iom_put( 'emp_ice', emp_ice )   ! emp over ice   (taking into account the snow blown away from the ice)
213
214      !                           ! vfxice = vfxbog + vfxbom + vfxsum + vfxsni + vfxopw + vfxdyn + vfxres + vfxlam + vfxpnd
215      CALL iom_put( 'vfxice'    , wfx_ice     )   ! mass flux from total ice growth/melt
216      CALL iom_put( 'vfxbog'    , wfx_bog     )   ! mass flux from bottom growth
217      CALL iom_put( 'vfxbom'    , wfx_bom     )   ! mass flux from bottom melt
218      CALL iom_put( 'vfxsum'    , wfx_sum     )   ! mass flux from surface melt
219      CALL iom_put( 'vfxlam'    , wfx_lam     )   ! mass flux from lateral melt
220      CALL iom_put( 'vfxsni'    , wfx_sni     )   ! mass flux from snow-ice formation
221      CALL iom_put( 'vfxopw'    , wfx_opw     )   ! mass flux from growth in open water
222      CALL iom_put( 'vfxdyn'    , wfx_dyn     )   ! mass flux from dynamics (ridging)
223      CALL iom_put( 'vfxres'    , wfx_res     )   ! mass flux from undiagnosed processes
224      CALL iom_put( 'vfxpnd'    , wfx_pnd     )   ! mass flux from melt ponds
225      CALL iom_put( 'vfxsub'    , wfx_ice_sub )   ! mass flux from ice sublimation (ice-atm.)
226      CALL iom_put( 'vfxsub_err', wfx_err_sub )   ! "excess" of sublimation sent to ocean     
227
228      IF ( iom_use( 'vfxthin' ) ) THEN   ! mass flux from ice growth in open water + thin ice (<20cm) => comparable to observations 
229         WHERE( hm_i(:,:) < 0.2 .AND. hm_i(:,:) > 0. ) ; z2d = wfx_bog
230         ELSEWHERE                                     ; z2d = 0._wp
231         END WHERE
232         CALL iom_put( 'vfxthin', wfx_opw + z2d )
233      ENDIF
234
235      !                            ! vfxsnw = vfxsnw_sni + vfxsnw_dyn + vfxsnw_sum
236      CALL iom_put( 'vfxsnw'     , wfx_snw     )   ! mass flux from total snow growth/melt
237      CALL iom_put( 'vfxsnw_sum' , wfx_snw_sum )   ! mass flux from snow melt at the surface
238      CALL iom_put( 'vfxsnw_sni' , wfx_snw_sni )   ! mass flux from snow melt during snow-ice formation
239      CALL iom_put( 'vfxsnw_dyn' , wfx_snw_dyn )   ! mass flux from dynamics (ridging)
240      CALL iom_put( 'vfxsnw_sub' , wfx_snw_sub )   ! mass flux from snow sublimation (ice-atm.)
241      CALL iom_put( 'vfxsnw_pre' , wfx_spr     )   ! snow precip
242
243      ! --- heat fluxes [W/m2] --- !
244      !                              ! qt_atm_oi - qt_oce_ai = hfxdhc - ( dihctrp + dshctrp )
245      IF( iom_use('qsr_oce'    ) )   CALL iom_put( 'qsr_oce'    , qsr_oce * ( 1._wp - at_i_b )                               )   !     solar flux at ocean surface
246      IF( iom_use('qns_oce'    ) )   CALL iom_put( 'qns_oce'    , qns_oce * ( 1._wp - at_i_b ) + qemp_oce                    )   ! non-solar flux at ocean surface
247      IF( iom_use('qsr_ice'    ) )   CALL iom_put( 'qsr_ice'    , SUM( qsr_ice * a_i_b, dim=3 )                              )   !     solar flux at ice surface
248      IF( iom_use('qns_ice'    ) )   CALL iom_put( 'qns_ice'    , SUM( qns_ice * a_i_b, dim=3 ) + qemp_ice                   )   ! non-solar flux at ice surface
249      IF( iom_use('qtr_ice_bot') )   CALL iom_put( 'qtr_ice_bot', SUM( qtr_ice_bot * a_i_b, dim=3 )                          )   !     solar flux transmitted thru ice
250      IF( iom_use('qtr_ice_top') )   CALL iom_put( 'qtr_ice_top', SUM( qtr_ice_top * a_i_b, dim=3 )                          )   !     solar flux transmitted thru ice surface
251      IF( iom_use('qt_oce'     ) )   CALL iom_put( 'qt_oce'     ,      ( qsr_oce + qns_oce ) * ( 1._wp - at_i_b ) + qemp_oce )
252      IF( iom_use('qt_ice'     ) )   CALL iom_put( 'qt_ice'     , SUM( ( qns_ice + qsr_ice ) * a_i_b, dim=3 )     + qemp_ice )
253      IF( iom_use('qt_oce_ai'  ) )   CALL iom_put( 'qt_oce_ai'  , qt_oce_ai * tmask(:,:,1)                                   )   ! total heat flux at the ocean   surface: interface oce-(ice+atm)
254      IF( iom_use('qt_atm_oi'  ) )   CALL iom_put( 'qt_atm_oi'  , qt_atm_oi * tmask(:,:,1)                                   )   ! total heat flux at the oce-ice surface: interface atm-(ice+oce)
255      IF( iom_use('qemp_oce'   ) )   CALL iom_put( 'qemp_oce'   , qemp_oce                                                   )   ! Downward Heat Flux from E-P over ocean
256      IF( iom_use('qemp_ice'   ) )   CALL iom_put( 'qemp_ice'   , qemp_ice                                                   )   ! Downward Heat Flux from E-P over ice
257
258      ! heat fluxes from ice transformations
259      !                            ! hfxdhc = hfxbog + hfxbom + hfxsum + hfxopw + hfxdif + hfxsnw - ( hfxthd + hfxdyn + hfxres + hfxsub + hfxspr )
260      CALL iom_put ('hfxbog'     , hfx_bog     )   ! heat flux used for ice bottom growth
261      CALL iom_put ('hfxbom'     , hfx_bom     )   ! heat flux used for ice bottom melt
262      CALL iom_put ('hfxsum'     , hfx_sum     )   ! heat flux used for ice surface melt
263      CALL iom_put ('hfxopw'     , hfx_opw     )   ! heat flux used for ice formation in open water
264      CALL iom_put ('hfxdif'     , hfx_dif     )   ! heat flux used for ice temperature change
265      CALL iom_put ('hfxsnw'     , hfx_snw     )   ! heat flux used for snow melt
266      CALL iom_put ('hfxerr'     , hfx_err_dif )   ! heat flux error after heat diffusion (included in qt_oce_ai)
267
268      ! heat fluxes associated with mass exchange (freeze/melt/precip...)
269      CALL iom_put ('hfxthd'     , hfx_thd     )   ! 
270      CALL iom_put ('hfxdyn'     , hfx_dyn     )   ! 
271      CALL iom_put ('hfxres'     , hfx_res     )   ! 
272      CALL iom_put ('hfxsub'     , hfx_sub     )   ! 
273      CALL iom_put ('hfxspr'     , hfx_spr     )   ! Heat flux from snow precip heat content
274
275      ! other heat fluxes
276      IF( iom_use('hfxsensib'  ) )   CALL iom_put( 'hfxsensib'  ,     -qsb_ice_bot * at_i_b         )   ! Sensible oceanic heat flux
277      IF( iom_use('hfxcndbot'  ) )   CALL iom_put( 'hfxcndbot'  , SUM( qcn_ice_bot * a_i_b, dim=3 ) )   ! Bottom conduction flux
278      IF( iom_use('hfxcndtop'  ) )   CALL iom_put( 'hfxcndtop'  , SUM( qcn_ice_top * a_i_b, dim=3 ) )   ! Surface conduction flux
279
280      ! controls
281      !---------
282#if ! defined key_agrif
283      IF( ln_icediachk .AND. .NOT. ln_bdy)   CALL ice_cons_final('iceupdate')                                       ! conservation
284#endif
285      IF( ln_icectl                      )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints
286      IF( sn_cfctl%l_prtctl              )   CALL ice_prt3D     ('iceupdate')                                       ! prints
287      IF( ln_timing                      )   CALL timing_stop   ('ice_update')                                      ! timing
288      !
289   END SUBROUTINE ice_update_flx
290
291
292   SUBROUTINE ice_update_tau( kt, pu_oce, pv_oce )
293      !!-------------------------------------------------------------------
294      !!                ***  ROUTINE ice_update_tau ***
295      !! 
296      !! ** Purpose : Update the ocean surface stresses due to the ice
297      !!         
298      !! ** Action  : * at each ice time step (every nn_fsbc time step):
299      !!                - compute the modulus of ice-ocean relative velocity
300      !!                  (*rho*Cd) at T-point (C-grid) or I-point (B-grid)
301      !!                      tmod_io = rhoco * | U_ice-U_oce |
302      !!                - update the modulus of stress at ocean surface
303      !!                      taum = (1-a) * taum + a * tmod_io * | U_ice-U_oce |
304      !!              * at each ocean time step (every kt):
305      !!                  compute linearized ice-ocean stresses as
306      !!                      Utau = tmod_io * | U_ice - pU_oce |
307      !!                using instantaneous current ocean velocity (usually before)
308      !!
309      !!    NB: - ice-ocean rotation angle no more allowed
310      !!        - here we make an approximation: taum is only computed every ice time step
311      !!          This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
312      !!          to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
313      !!
314      !! ** Outputs : - utau, vtau   : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
315      !!              - taum         : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
316      !!---------------------------------------------------------------------
317      INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index
318      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents
319      !
320      INTEGER  ::   ji, jj   ! dummy loop indices
321      REAL(wp) ::   zat_u, zutau_ice, zu_t, zmodt   ! local scalar
322      REAL(wp) ::   zat_v, zvtau_ice, zv_t, zrhoco  !   -      -
323      !!---------------------------------------------------------------------
324      IF( ln_timing )   CALL timing_start('ice_update_tau')
325
326      IF( kt == nit000 .AND. lwp ) THEN
327         WRITE(numout,*)
328         WRITE(numout,*)'ice_update_tau: update stress at the ice-ocean interface'
329         WRITE(numout,*)'~~~~~~~~~~~~~~'
330      ENDIF
331
332      zrhoco = rau0 * rn_cio
333      !
334      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
335         DO_2D_00_00
336            !                                               ! 2*(U_ice-U_oce) at T-point
337            zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)   
338            zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) 
339            !                                              ! |U_ice-U_oce|^2
340            zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  )
341            !                                               ! update the ocean stress modulus
342            taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * zrhoco * zmodt
343            tmod_io(ji,jj) = zrhoco * SQRT( zmodt )          ! rhoco * |U_ice-U_oce| at T-point
344         END_2D
345         CALL lbc_lnk_multi( 'iceupdate', taum, 'T', 1., tmod_io, 'T', 1. )
346         !
347         utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step
348         vtau_oce(:,:) = vtau(:,:)
349         !
350      ENDIF
351      !
352      !                                      !==  every ocean time-step  ==!
353      !
354      DO_2D_00_00
355         ! ice area at u and v-points
356         zat_u  = ( at_i(ji,jj) * tmask(ji,jj,1) + at_i (ji+1,jj    ) * tmask(ji+1,jj  ,1) )  &
357            &     / MAX( 1.0_wp , tmask(ji,jj,1) + tmask(ji+1,jj  ,1) )
358         zat_v  = ( at_i(ji,jj) * tmask(ji,jj,1) + at_i (ji  ,jj+1  ) * tmask(ji  ,jj+1,1) )  &
359            &     / MAX( 1.0_wp , tmask(ji,jj,1) + tmask(ji  ,jj+1,1) )
360         !                                                   ! linearized quadratic drag formulation
361         zutau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
362         zvtau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
363         !                                                   ! stresses at the ocean surface
364         utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
365         vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
366      END_2D
367      CALL lbc_lnk_multi( 'iceupdate', utau, 'U', -1., vtau, 'V', -1. )   ! lateral boundary condition
368      !
369      IF( ln_timing )   CALL timing_stop('ice_update_tau')
370      ! 
371   END SUBROUTINE ice_update_tau
372
373
374   SUBROUTINE ice_update_init
375      !!-------------------------------------------------------------------
376      !!                  ***  ROUTINE ice_update_init  ***
377      !!             
378      !! ** Purpose :   allocate ice-ocean stress fields and read restarts
379      !!                containing the snow & ice mass
380      !!
381      !!-------------------------------------------------------------------
382      INTEGER  ::   ji, jj, jk               ! dummy loop indices
383      REAL(wp) ::   zcoefu, zcoefv, zcoeff   ! local scalar
384      !!-------------------------------------------------------------------
385      !
386      IF(lwp) WRITE(numout,*)
387      IF(lwp) WRITE(numout,*) 'ice_update_init: ice-ocean stress init'
388      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
389      !
390      !                                      ! allocate ice_update array
391      IF( ice_update_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_update_init : unable to allocate standard arrays' )
392      !
393      CALL update_rst( 'READ' )  !* read or initialize all required files
394      !
395   END SUBROUTINE ice_update_init
396
397
398   SUBROUTINE update_rst( cdrw, kt )
399      !!---------------------------------------------------------------------
400      !!                   ***  ROUTINE rhg_evp_rst  ***
401      !!                     
402      !! ** Purpose :   Read or write RHG file in restart file
403      !!
404      !! ** Method  :   use of IOM library
405      !!----------------------------------------------------------------------
406      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! 'READ'/'WRITE' flag
407      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
408      !
409      INTEGER  ::   iter   ! local integer
410      INTEGER  ::   id1    ! local integer
411      !!----------------------------------------------------------------------
412      !
413      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
414         !                                   ! ---------------
415         IF( ln_rstart ) THEN                   !* Read the restart file
416            !
417            id1 = iom_varid( numrir, 'snwice_mass' , ldstop = .FALSE. )
418            !
419            IF( id1 > 0 ) THEN                       ! fields exist
420               CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass'  , snwice_mass   )
421               CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b )
422            ELSE                                     ! start from rest
423               IF(lwp) WRITE(numout,*) '   ==>>   previous run without snow-ice mass output then set it'
424               snwice_mass  (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) )
425               snwice_mass_b(:,:) = snwice_mass(:,:)
426            ENDIF
427         ELSE                                   !* Start from rest
428            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set the snow-ice mass'
429            snwice_mass  (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) )
430            snwice_mass_b(:,:) = snwice_mass(:,:)
431         ENDIF
432         !
433      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
434         !                                   ! -------------------
435         IF(lwp) WRITE(numout,*) '---- update-rst ----'
436         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
437         !
438         CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass'  , snwice_mass   )
439         CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b )
440         !
441      ENDIF
442      !
443   END SUBROUTINE update_rst
444
445#else
446   !!----------------------------------------------------------------------
447   !!   Default option         Dummy module           NO SI3 sea-ice model
448   !!----------------------------------------------------------------------
449#endif 
450
451   !!======================================================================
452END MODULE iceupdate
Note: See TracBrowser for help on using the repository browser.