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_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90 @ 8879

Last change on this file since 8879 was 8879, checked in by frrh, 6 years ago

Merge in http://fcm3/projects/NEMO.xm/log/branches/UKMO/dev_r8183_ICEMODEL_svn_removed
revisions 8738:8847 inclusive.

File size: 25.7 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) + wfx_pnd(ji,jj)
167
168            ! add the snow melt water to snow mass flux to the ocean
169            wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj)
170
171            ! mass flux at the ocean/ice interface
172            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
173            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)
174
175
176            ! Salt flux at the ocean surface     
177            !------------------------------------------
178            sfx(ji,jj) = sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj)   &
179               &       + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj)
180           
181            ! Mass of snow and ice per unit area   
182            !----------------------------------------
183            snwice_mass_b(ji,jj) = snwice_mass(ji,jj)       ! save mass from the previous ice time step
184            !                                               ! new mass per unit area
185            snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj)  ) 
186            !                                               ! time evolution of snow+ice mass
187            snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_rdtice
188           
189         END DO
190      END DO
191
192      ! Storing the transmitted variables
193      !----------------------------------
194      fr_i  (:,:)   = at_i(:,:)             ! Sea-ice fraction           
195      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                     
196
197      ! Snow/ice albedo (only if sent to coupler, useless in forced mode)
198      !------------------------------------------------------------------
199      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
200      !
201      alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
202      !
203      IF( lrst_ice ) THEN                       !* write snwice_mass fields in the restart file
204         CALL update_rst( 'WRITE', kt )
205      ENDIF
206      !
207      ! output all fluxes
208      !------------------
209      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * ( 1._wp - at_i_b(:,:) ) )                      !     solar flux at ocean surface
210      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
211      IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface
212      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
213      IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice
214      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) )
215      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   &
216         &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) )
217!!gm I don't understand the variable below.... why not multiplied by a_i_b or (1-a_i_b) ???
218      IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) ) 
219      IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) ) 
220      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)
221      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)
222
223      CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation [m/day]
224
225      CALL iom_put( "sfxbog"      , sfx_bog             )        ! salt flux from bottom growth
226      CALL iom_put( "sfxbom"      , sfx_bom             )        ! salt flux from bottom melting
227      CALL iom_put( "sfxsum"      , sfx_sum             )        ! salt flux from surface melting
228      CALL iom_put( "sfxlam"      , sfx_lam             )        ! salt flux from lateral melting
229      CALL iom_put( "sfxsni"      , sfx_sni             )        ! salt flux from snow ice formation
230      CALL iom_put( "sfxopw"      , sfx_opw             )        ! salt flux from open water formation
231      CALL iom_put( "sfxdyn"      , sfx_dyn             )        ! salt flux from ridging rafting
232      CALL iom_put( "sfxres"      , sfx_res             )        ! salt flux from corrections (resultant)
233      CALL iom_put( "sfxbri"      , sfx_bri             )        ! salt flux from brines
234      CALL iom_put( "sfxsub"      , sfx_sub             )        ! salt flux from sublimation
235      CALL iom_put( "sfx"         , sfx                 )        ! total salt flux
236
237      CALL iom_put( "vfxres"     , wfx_res              )        ! prod./melting due to corrections
238      CALL iom_put( "vfxopw"     , wfx_opw              )        ! lateral thermodynamic ice production
239      CALL iom_put( "vfxsni"     , wfx_sni              )        ! snowice ice production
240      CALL iom_put( "vfxbog"     , wfx_bog              )        ! bottom thermodynamic ice production
241      CALL iom_put( "vfxdyn"     , wfx_dyn              )        ! dynamic ice production (rid/raft)
242      CALL iom_put( "vfxsum"     , wfx_sum              )        ! surface melt
243      CALL iom_put( "vfxbom"     , wfx_bom              )        ! bottom melt
244      CALL iom_put( "vfxlam"     , wfx_lam              )        ! lateral melt
245      CALL iom_put( "vfxice"     , wfx_ice              )        ! total ice growth/melt
246      IF ( ln_pnd_fwb )   CALL iom_put( "vfxpnd", wfx_pnd   )        ! melt pond water flux
247
248      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations 
249         WHERE( hm_i(:,:) < 0.2 .AND. hm_i(:,:) > 0. ) ; z2d = wfx_bog
250         ELSEWHERE                                     ; z2d = 0._wp
251         END WHERE
252         CALL iom_put( "vfxthin", ( wfx_opw + z2d )     )
253      ENDIF
254
255      CALL iom_put( "vfxspr"     , wfx_spr              )        ! precip (snow)
256      CALL iom_put( "vfxsnw"     , wfx_snw              )        ! total snw growth/melt
257      CALL iom_put( "vfxsub"     , wfx_sub              )        ! sublimation (snow/ice)
258      CALL iom_put( "vfxsub_err" , wfx_err_sub          )        ! "excess" of sublimation sent to ocean     
259 
260      CALL iom_put ('hfxthd'     , hfx_thd(:,:)         )   
261      CALL iom_put ('hfxdyn'     , hfx_dyn(:,:)         )   
262      CALL iom_put ('hfxres'     , hfx_res(:,:)         )   
263      CALL iom_put ('hfxout'     , hfx_out(:,:)         )   
264      CALL iom_put ('hfxin'      , hfx_in(:,:)          )   
265      CALL iom_put ('hfxsnw'     , hfx_snw(:,:)         )   
266      CALL iom_put ('hfxsub'     , hfx_sub(:,:)         )   
267      CALL iom_put ('hfxerr'     , hfx_err_dif(:,:)     )   
268      CALL iom_put ('hfxerr_rem' , hfx_err_rem(:,:)     )   
269     
270      CALL iom_put ('hfxsum'     , hfx_sum(:,:)         )   
271      CALL iom_put ('hfxbom'     , hfx_bom(:,:)         )   
272      CALL iom_put ('hfxbog'     , hfx_bog(:,:)         )   
273      CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   
274      CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   
275      CALL iom_put ('hfxtur'     , fhtur(:,:) * at_i_b(:,:) ) ! turbulent heat flux at ice base
276      CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice
277      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip
278      !
279      ! controls
280      !---------
281      IF( ln_icediachk .AND. .NOT. ln_bdy)   CALL ice_cons_final('iceupdate')                                       ! conservation
282      IF( ln_icectl                      )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints
283      IF( ln_ctl                         )   CALL ice_prt3D     ('iceupdate')                                       ! prints
284      IF( nn_timing == 1                 )   CALL timing_stop   ('ice_update')                                      ! timing
285      !
286   END SUBROUTINE ice_update_flx
287
288
289   SUBROUTINE ice_update_tau( kt, pu_oce, pv_oce )
290      !!-------------------------------------------------------------------
291      !!                ***  ROUTINE ice_update_tau ***
292      !! 
293      !! ** Purpose : Update the ocean surface stresses due to the ice
294      !!         
295      !! ** Action  : * at each ice time step (every nn_fsbc time step):
296      !!                - compute the modulus of ice-ocean relative velocity
297      !!                  (*rho*Cd) at T-point (C-grid) or I-point (B-grid)
298      !!                      tmod_io = rhoco * | U_ice-U_oce |
299      !!                - update the modulus of stress at ocean surface
300      !!                      taum = (1-a) * taum + a * tmod_io * | U_ice-U_oce |
301      !!              * at each ocean time step (every kt):
302      !!                  compute linearized ice-ocean stresses as
303      !!                      Utau = tmod_io * | U_ice - pU_oce |
304      !!                using instantaneous current ocean velocity (usually before)
305      !!
306      !!    NB: - ice-ocean rotation angle no more allowed
307      !!        - here we make an approximation: taum is only computed every ice time step
308      !!          This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
309      !!          to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
310      !!
311      !! ** Outputs : - utau, vtau   : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
312      !!              - taum         : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
313      !!---------------------------------------------------------------------
314      INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index
315      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents
316      !
317      INTEGER  ::   ji, jj   ! dummy loop indices
318      REAL(wp) ::   zat_u, zutau_ice, zu_t, zmodt   ! local scalar
319      REAL(wp) ::   zat_v, zvtau_ice, zv_t, zrhoco  !   -      -
320      !!---------------------------------------------------------------------
321
322      IF( nn_timing == 1 )  CALL timing_start('ice_update')
323
324      IF( kt == nit000 .AND. lwp ) THEN
325         WRITE(numout,*)
326         WRITE(numout,*)'ice_update_tau: update stress at the ice-ocean interface'
327         WRITE(numout,*)'~~~~~~~~~~~~~~'
328      ENDIF
329
330      zrhoco = rau0 * rn_cio
331      !
332      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
333         DO jj = 2, jpjm1                             !* update the modulus of stress at ocean surface (T-point)
334            DO ji = fs_2, fs_jpim1
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 DO
344         END DO
345         CALL lbc_lnk_multi( 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 jj = 2, jpjm1                                !* update the stress WITHOUT an ice-ocean rotation angle
355         DO ji = fs_2, fs_jpim1   ! Vect. Opt.
356            zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp   ! ice area at u and V-points
357            zat_v  = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5_wp
358            !                                                   ! linearized quadratic drag formulation
359            zutau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
360            zvtau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
361            !                                                   ! stresses at the ocean surface
362            utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
363            vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
364         END DO
365      END DO
366      CALL lbc_lnk_multi( utau, 'U', -1., vtau, 'V', -1. )   ! lateral boundary condition
367      !
368      IF( nn_timing == 1 )  CALL timing_stop('ice_update')
369     
370   END SUBROUTINE ice_update_tau
371
372
373   SUBROUTINE ice_update_init
374      !!-------------------------------------------------------------------
375      !!                  ***  ROUTINE ice_update_init  ***
376      !!             
377      !! ** Purpose :   ???
378      !!
379      !!-------------------------------------------------------------------
380      INTEGER  ::   ji, jj, jk               ! dummy loop indices
381      REAL(wp) ::   zcoefu, zcoefv, zcoeff   ! local scalar
382      !!-------------------------------------------------------------------
383      !
384      IF(lwp) WRITE(numout,*)
385      IF(lwp) WRITE(numout,*) 'ice_update_init: ???? '
386      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
387
388      !                                      ! allocate ice_update array
389      IF( ice_update_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_update_init : unable to allocate standard arrays' )
390      !
391      CALL update_rst( 'READ' )  !* read or initialize all required files
392      !
393   END SUBROUTINE ice_update_init
394
395   SUBROUTINE update_rst( cdrw, kt )
396      !!---------------------------------------------------------------------
397      !!                   ***  ROUTINE rhg_evp_rst  ***
398      !!                     
399      !! ** Purpose :   Read or write RHG file in restart file
400      !!
401      !! ** Method  :   use of IOM library
402      !!----------------------------------------------------------------------
403      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
404      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
405      !
406      INTEGER  ::   iter   ! local integer
407      INTEGER  ::   id1    ! local integer
408      !!----------------------------------------------------------------------
409      !
410      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
411         !                                   ! ---------------
412         IF( ln_rstart ) THEN                   !* Read the restart file
413            !
414            id1 = iom_varid( numrir, 'snwice_mass' , ldstop = .FALSE. )
415            !
416            IF( id1 > 0 ) THEN                       ! fields exist
417               CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass'  , snwice_mass   )
418               CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b )
419            ELSE                                     ! start from rest
420               IF(lwp) WRITE(numout,*) '   ==>>   previous run without snow-ice mass output then set it'
421               snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
422               snwice_mass_b(:,:) = snwice_mass(:,:)
423            ENDIF
424         ELSE                                   !* Start from rest
425            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set the snow-ice mass'
426            snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
427            snwice_mass_b(:,:) = snwice_mass(:,:)
428         ENDIF
429         !
430      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
431         !                                   ! -------------------
432         IF(lwp) WRITE(numout,*) '---- update-rst ----'
433         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
434         !
435         CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass'  , snwice_mass   )
436         CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b )
437         !
438      ENDIF
439      !
440   END SUBROUTINE update_rst
441
442#else
443   !!----------------------------------------------------------------------
444   !!   Default option         Dummy module           NO ESIM sea-ice model
445   !!----------------------------------------------------------------------
446#endif 
447
448   !!======================================================================
449END MODULE iceupdate
Note: See TracBrowser for help on using the repository browser.