source: NEMO/trunk/src/ICE/iceupdate.F90 @ 12489

Last change on this file since 12489 was 12489, checked in by davestorkey, 7 months ago

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp → rn_atfp (namelist parameter used everywhere) and rau0 → rho0.

This version gives bit-comparable results to the previous version of the trunk.

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