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 branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90 @ 8738

Last change on this file since 8738 was 8738, checked in by dancopsey, 6 years ago

Merged in main ICEMODEL branch (branches/2017/dev_r8183_ICEMODEL) up to revision 8588

File size: 25.8 KB
Line 
1MODULE iceupdate
2   !!======================================================================
3   !!                       ***  MODULE iceupdate   ***
4   !!  Sea-ice :   computation of the flux at the sea ice/ocean interface
5   !!======================================================================
6   !! History :   -   ! 2006-07 (M. Vancoppelle)  LIM3 original code
7   !!            3.0  ! 2008-03 (C. Tallandier)  surface module
8   !!             -   ! 2008-04 (C. Tallandier)  split in 2 + new ice-ocean coupling
9   !!            3.3  ! 2010-05 (G. Madec) decrease ocean & ice reference salinities in the Baltic sea
10   !!                 !                  + simplification of the ice-ocean stress calculation
11   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
12   !!             -   ! 2012    (D. Iovino) salt flux change
13   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
14   !!            3.5  ! 2012-10 (A. Coward, G. Madec) salt fluxes ; ice+snow mass
15   !!----------------------------------------------------------------------
16#if defined key_lim3
17   !!----------------------------------------------------------------------
18   !!   'key_lim3'                                       ESIM sea-ice model
19   !!----------------------------------------------------------------------
20   !!   ice_update_alloc : allocate the iceupdate arrays
21   !!   ice_update_init  : initialisation
22   !!   ice_update_flx   : updates mass, heat and salt fluxes at the ocean surface
23   !!   ice_update_tau   : update i- and j-stresses, and its modulus at the ocean surface
24   !!----------------------------------------------------------------------
25   USE oce     , ONLY : sshn, sshb
26   USE phycst         ! physical constants
27   USE dom_oce        ! ocean domain
28   USE ice            ! sea-ice: variables
29!!gm  It should be probably better to pass these variable in argument of the routine,
30!!gm  rather than having this long list in USE. This will also highlight what is updated, and what is just use.
31   USE sbc_ice
32   USE sbc_oce
33!!gm end
34   USE sbccpl         ! Surface boundary condition: coupled interface
35   USE icealb         ! sea-ice: albedo parameters
36   USE traqsr         ! add penetration of solar flux in the calculation of heat budget
37   USE icectl         ! sea-ice: control prints
38   USE bdy_oce , ONLY : ln_bdy
39   !
40   USE in_out_manager ! I/O manager
41   USE iom            ! I/O manager library
42   USE lib_mpp        ! MPP library
43   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
44   USE lbclnk         ! lateral boundary conditions (or mpp links)
45   USE timing         ! Timing
46
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC   ice_update_init   ! called by ice_init
51   PUBLIC   ice_update_flx    ! called by ice_stp
52   PUBLIC   ice_update_tau    ! called by ice_stp
53
54   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress     [N/m2]
55   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmod_io              ! modulus of the ice-ocean velocity   [m/s]
56
57   !! * Substitutions
58#  include "vectopt_loop_substitute.h90"
59   !!----------------------------------------------------------------------
60   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
61   !! $Id: iceupdate.F90 8411 2017-08-07 16:09:12Z clem $
62   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   INTEGER FUNCTION ice_update_alloc()
67      !!-------------------------------------------------------------------
68      !!             ***  ROUTINE ice_update_alloc ***
69      !!-------------------------------------------------------------------
70      ALLOCATE( utau_oce(jpi,jpj), vtau_oce(jpi,jpj), tmod_io(jpi,jpj), STAT=ice_update_alloc )
71      !
72      IF( lk_mpp                )   CALL mpp_sum( ice_update_alloc )
73      IF( ice_update_alloc /= 0 )   CALL ctl_warn('ice_update_alloc: failed to allocate arrays')
74      !
75   END FUNCTION ice_update_alloc
76
77
78   SUBROUTINE ice_update_flx( kt )
79      !!-------------------------------------------------------------------
80      !!                ***  ROUTINE ice_update_flx ***
81      !! 
82      !! ** Purpose :   Update the surface ocean boundary condition for heat
83      !!                salt and mass over areas where sea-ice is non-zero
84      !!         
85      !! ** Action  : - computes the heat and freshwater/salt fluxes
86      !!                at the ice-ocean interface.
87      !!              - Update the ocean sbc
88      !!     
89      !! ** Outputs : - qsr     : sea heat flux:     solar
90      !!              - qns     : sea heat flux: non solar
91      !!              - emp     : freshwater budget: volume flux
92      !!              - sfx     : salt flux
93      !!              - fr_i    : ice fraction
94      !!              - tn_ice  : sea-ice surface temperature
95      !!              - alb_ice : sea-ice albedo (recomputed only for coupled mode)
96      !!
97      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
98      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
99      !!              These refs are now obsolete since everything has been revised
100      !!              The ref should be Rousset et al., 2015
101      !!---------------------------------------------------------------------
102      INTEGER, INTENT(in) ::   kt   ! number of iteration
103      !
104      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
105      REAL(wp) ::   zqmass           ! Heat flux associated with mass exchange ice->ocean (W.m-2)
106      REAL(wp) ::   zqsr             ! New solar flux received by the ocean
107      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d                  ! 2D workspace
108      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_cs, zalb_os     ! 3D workspace
109      !!---------------------------------------------------------------------
110      IF( nn_timing == 1 )  CALL timing_start('ice_update')
111
112      IF( kt == nit000 .AND. lwp ) THEN
113         WRITE(numout,*)
114         WRITE(numout,*)'ice_update_flx: update fluxes (mass, salt and heat) at the ice-ocean interface'
115         WRITE(numout,*)'~~~~~~~~~~~~~~'
116      ENDIF
117
118      ! --- case we bypass ice thermodynamics --- !
119      IF( .NOT. ln_icethd ) THEN   ! we suppose ice is impermeable => ocean is isolated from atmosphere
120         hfx_in   (:,:)   = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
121         hfx_out  (:,:)   = ( 1._wp - at_i_b(:,:) ) *   qns_oce(:,:)                  + qemp_oce(:,:)
122         ftr_ice  (:,:,:) = 0._wp
123         emp_ice  (:,:)   = 0._wp
124         qemp_ice (:,:)   = 0._wp
125         qevap_ice(:,:,:) = 0._wp
126      ENDIF
127     
128      DO jj = 1, jpj
129         DO ji = 1, jpi
130
131            ! Solar heat flux reaching the ocean = zqsr (W.m-2)
132            !---------------------------------------------------
133            zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - ftr_ice(ji,jj,:) ) )
134
135            ! Total heat flux reaching the ocean = hfx_out (W.m-2)
136            !---------------------------------------------------
137            zqmass         = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC)
138            hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr
139
140            ! Add the residual from heat diffusion equation and sublimation (W.m-2)
141            !----------------------------------------------------------------------
142            hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj) +   &
143               &           ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) )
144
145            ! New qsr and qns used to compute the oceanic heat flux at the next time step
146            !----------------------------------------------------------------------------
147            qsr(ji,jj) = zqsr                                     
148            qns(ji,jj) = hfx_out(ji,jj) - zqsr             
149
150            ! Mass flux at the atm. surface       
151            !-----------------------------------
152            wfx_sub(ji,jj) = wfx_snw_sub(ji,jj) + wfx_ice_sub(ji,jj)
153
154            ! Mass flux at the ocean surface     
155            !------------------------------------
156            !  case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED)
157            !  -------------------------------------------------------------------------------------
158            !  The idea of this approach is that the system that we consider is the ICE-OCEAN system
159            !  Thus  FW  flux  =  External ( E-P+snow melt)
160            !       Salt flux  =  Exchanges in the ice-ocean system then converted into FW
161            !                     Associated to Ice formation AND Ice melting
162            !                     Even if i see Ice melting as a FW and SALT flux
163            !       
164            ! mass flux from ice/ocean
165            wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   &
166               &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) 
167
168            IF ( ln_pnd_fw )   wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj)
169
170            ! add the snow melt water to snow mass flux to the ocean
171            wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj)
172
173            ! mass flux at the ocean/ice interface
174            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
175            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)
176
177
178            ! Salt flux at the ocean surface     
179            !------------------------------------------
180            sfx(ji,jj) = sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj)   &
181               &       + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj)
182           
183            ! Mass of snow and ice per unit area   
184            !----------------------------------------
185            snwice_mass_b(ji,jj) = snwice_mass(ji,jj)       ! save mass from the previous ice time step
186            !                                               ! new mass per unit area
187            snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj)  ) 
188            !                                               ! time evolution of snow+ice mass
189            snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_rdtice
190           
191         END DO
192      END DO
193
194      ! Storing the transmitted variables
195      !----------------------------------
196      fr_i  (:,:)   = at_i(:,:)             ! Sea-ice fraction           
197      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                     
198
199      ! Snow/ice albedo (only if sent to coupler, useless in forced mode)
200      !------------------------------------------------------------------
201      CALL ice_alb( t_su, h_i, h_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
202      !
203      alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
204      !
205      IF( lrst_ice ) THEN                       !* write snwice_mass fields in the restart file
206         CALL update_rst( 'WRITE', kt )
207      ENDIF
208      !
209      ! output all fluxes
210      !------------------
211      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * ( 1._wp - at_i_b(:,:) ) )                      !     solar flux at ocean surface
212      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
213      IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface
214      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
215      IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice
216      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) )
217      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   &
218         &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) )
219!!gm I don't understand the variable below.... why not multiplied by a_i_b or (1-a_i_b) ???
220      IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) ) 
221      IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) ) 
222      IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce (:,:) )   ! emp over ocean (taking into account the snow blown away from the ice)
223      IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice (:,:) )   ! emp over ice   (taking into account the snow blown away from the ice)
224
225      CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation [m/day]
226
227      CALL iom_put( "sfxbog"      , sfx_bog             )        ! salt flux from bottom growth
228      CALL iom_put( "sfxbom"      , sfx_bom             )        ! salt flux from bottom melting
229      CALL iom_put( "sfxsum"      , sfx_sum             )        ! salt flux from surface melting
230      CALL iom_put( "sfxlam"      , sfx_lam             )        ! salt flux from lateral melting
231      CALL iom_put( "sfxsni"      , sfx_sni             )        ! salt flux from snow ice formation
232      CALL iom_put( "sfxopw"      , sfx_opw             )        ! salt flux from open water formation
233      CALL iom_put( "sfxdyn"      , sfx_dyn             )        ! salt flux from ridging rafting
234      CALL iom_put( "sfxres"      , sfx_res             )        ! salt flux from corrections (resultant)
235      CALL iom_put( "sfxbri"      , sfx_bri             )        ! salt flux from brines
236      CALL iom_put( "sfxsub"      , sfx_sub             )        ! salt flux from sublimation
237      CALL iom_put( "sfx"         , sfx                 )        ! total salt flux
238
239      CALL iom_put( "vfxres"     , wfx_res              )        ! prod./melting due to corrections
240      CALL iom_put( "vfxopw"     , wfx_opw              )        ! lateral thermodynamic ice production
241      CALL iom_put( "vfxsni"     , wfx_sni              )        ! snowice ice production
242      CALL iom_put( "vfxbog"     , wfx_bog              )        ! bottom thermodynamic ice production
243      CALL iom_put( "vfxdyn"     , wfx_dyn              )        ! dynamic ice production (rid/raft)
244      CALL iom_put( "vfxsum"     , wfx_sum              )        ! surface melt
245      CALL iom_put( "vfxbom"     , wfx_bom              )        ! bottom melt
246      CALL iom_put( "vfxlam"     , wfx_lam              )        ! lateral melt
247      CALL iom_put( "vfxice"     , wfx_ice              )        ! total ice growth/melt
248      IF ( ln_pnd )   CALL iom_put( "vfxpnd", wfx_pnd   )        ! melt pond water flux
249
250      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations 
251         WHERE( hm_i(:,:) < 0.2 .AND. hm_i(:,:) > 0. ) ; z2d = wfx_bog
252         ELSEWHERE                                     ; z2d = 0._wp
253         END WHERE
254         CALL iom_put( "vfxthin", ( wfx_opw + z2d )     )
255      ENDIF
256
257      CALL iom_put( "vfxspr"     , wfx_spr              )        ! precip (snow)
258      CALL iom_put( "vfxsnw"     , wfx_snw              )        ! total snw growth/melt
259      CALL iom_put( "vfxsub"     , wfx_sub              )        ! sublimation (snow/ice)
260      CALL iom_put( "vfxsub_err" , wfx_err_sub          )        ! "excess" of sublimation sent to ocean     
261 
262      CALL iom_put ('hfxthd'     , hfx_thd(:,:)         )   
263      CALL iom_put ('hfxdyn'     , hfx_dyn(:,:)         )   
264      CALL iom_put ('hfxres'     , hfx_res(:,:)         )   
265      CALL iom_put ('hfxout'     , hfx_out(:,:)         )   
266      CALL iom_put ('hfxin'      , hfx_in(:,:)          )   
267      CALL iom_put ('hfxsnw'     , hfx_snw(:,:)         )   
268      CALL iom_put ('hfxsub'     , hfx_sub(:,:)         )   
269      CALL iom_put ('hfxerr'     , hfx_err_dif(:,:)     )   
270      CALL iom_put ('hfxerr_rem' , hfx_err_rem(:,:)     )   
271     
272      CALL iom_put ('hfxsum'     , hfx_sum(:,:)         )   
273      CALL iom_put ('hfxbom'     , hfx_bom(:,:)         )   
274      CALL iom_put ('hfxbog'     , hfx_bog(:,:)         )   
275      CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   
276      CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   
277      CALL iom_put ('hfxtur'     , fhtur(:,:) * at_i_b(:,:) ) ! turbulent heat flux at ice base
278      CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice
279      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip
280      !
281      ! controls
282      !---------
283      IF( ln_icediachk .AND. .NOT. ln_bdy)   CALL ice_cons_final('iceupdate')                                       ! conservation
284      IF( ln_icectl                      )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints
285      IF( ln_ctl                         )   CALL ice_prt3D     ('iceupdate')                                       ! prints
286      IF( nn_timing == 1                 )   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
324      IF( nn_timing == 1 )  CALL timing_start('ice_update')
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 jj = 2, jpjm1                             !* update the modulus of stress at ocean surface (T-point)
336            DO ji = fs_2, fs_jpim1
337               !                                               ! 2*(U_ice-U_oce) at T-point
338               zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)   
339               zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) 
340               !                                              ! |U_ice-U_oce|^2
341               zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  )
342               !                                               ! update the ocean stress modulus
343               taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * zrhoco * zmodt
344               tmod_io(ji,jj) = zrhoco * SQRT( zmodt )          ! rhoco * |U_ice-U_oce| at T-point
345            END DO
346         END DO
347         CALL lbc_lnk_multi( taum, 'T', 1., tmod_io, 'T', 1. )
348         !
349         utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step
350         vtau_oce(:,:) = vtau(:,:)
351         !
352      ENDIF
353      !
354      !                                      !==  every ocean time-step  ==!
355      !
356      DO jj = 2, jpjm1                                !* update the stress WITHOUT an ice-ocean rotation angle
357         DO ji = fs_2, fs_jpim1   ! Vect. Opt.
358            zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp   ! ice area at u and V-points
359            zat_v  = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5_wp
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 DO
367      END DO
368      CALL lbc_lnk_multi( utau, 'U', -1., vtau, 'V', -1. )   ! lateral boundary condition
369      !
370      IF( nn_timing == 1 )  CALL timing_stop('ice_update')
371     
372   END SUBROUTINE ice_update_tau
373
374
375   SUBROUTINE ice_update_init
376      !!-------------------------------------------------------------------
377      !!                  ***  ROUTINE ice_update_init  ***
378      !!             
379      !! ** Purpose :   ???
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: ???? '
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   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) * ( rhosn * vt_s(:,:) + rhoic * 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) * ( rhosn * vt_s(:,:) + rhoic * 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 ESIM sea-ice model
447   !!----------------------------------------------------------------------
448#endif 
449
450   !!======================================================================
451END MODULE iceupdate
Note: See TracBrowser for help on using the repository browser.