source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90 @ 8486

Last change on this file since 8486 was 8486, checked in by clem, 3 years ago

changes in style - part1 - (now the code looks better txs to Gurvan's comments)

File size: 22.5 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'                                    LIM 3.0 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 par_oce        ! ocean parameters
26   USE oce     , ONLY : sshn, sshb
27   USE phycst         ! physical constants
28   USE dom_oce        ! ocean domain
29   USE ice            ! sea-ice: variables
30!!gm  It should be probably better to pass these variable in argument of the routine,
31!!gm  rather than having this long list in USE. This will also highlight what is updated, and what is just use.
32   USE sbc_ice , ONLY : emp_oce, qns_oce, qsr_oce , qemp_oce ,                             &
33      &                 emp_ice, qsr_ice, qemp_ice, qevap_ice, alb_ice, tn_ice, cldf_ice,  &
34      &                 snwice_mass, snwice_mass_b, snwice_fmass
35   USE sbc_oce , ONLY : nn_fsbc, ln_ice_embd, sfx, fr_i, qsr_tot, qns, qsr, fmmflx, emp, taum, utau, vtau
36!!gm end
37   USE sbccpl         ! Surface boundary condition: coupled interface
38   USE icealb         ! albedo parameters
39   USE traqsr         ! add penetration of solar flux in the calculation of heat budget
40   USE domvvl         ! Variable volume
41   USE icectl         ! ???
42   USE bdy_oce , ONLY : ln_bdy
43   !
44   USE in_out_manager ! I/O manager
45   USE iom            ! xIO server
46   USE lbclnk         ! ocean lateral boundary condition - MPP exchanges
47   USE lib_mpp        ! MPP library
48   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
49   USE timing         ! Timing
50
51   IMPLICIT NONE
52   PRIVATE
53
54   PUBLIC   ice_update_init   ! called by ice_init
55   PUBLIC   ice_update_flx    ! called by ice_stp
56   PUBLIC   ice_update_tau    ! called by ice_stp
57
58   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress     [N/m2]
59   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmod_io              ! modulus of the ice-ocean velocity   [m/s]
60   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   soce_0  , sice_0     ! cst SSS and ice salinity (levitating sea-ice)
61
62   !! * Substitutions
63#  include "vectopt_loop_substitute.h90"
64   !!----------------------------------------------------------------------
65   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
66   !! $Id: iceupdate.F90 8411 2017-08-07 16:09:12Z clem $
67   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
68   !!----------------------------------------------------------------------
69CONTAINS
70
71   INTEGER FUNCTION ice_update_alloc()
72      !!-------------------------------------------------------------------
73      !!             ***  ROUTINE ice_update_alloc ***
74      !!-------------------------------------------------------------------
75      ALLOCATE( soce_0(jpi,jpj) , utau_oce(jpi,jpj) ,                       &
76         &      sice_0(jpi,jpj) , vtau_oce(jpi,jpj) , tmod_io(jpi,jpj), STAT=ice_update_alloc)
77         !
78      IF( lk_mpp                )   CALL mpp_sum( ice_update_alloc )
79      IF( ice_update_alloc /= 0 )   CALL ctl_warn('ice_update_alloc: failed to allocate arrays')
80   END FUNCTION ice_update_alloc
81
82
83   SUBROUTINE ice_update_flx( kt )
84      !!-------------------------------------------------------------------
85      !!                ***  ROUTINE ice_update_flx ***
86      !! 
87      !! ** Purpose :   Update the surface ocean boundary condition for heat
88      !!              salt and mass over areas where sea-ice is non-zero
89      !!         
90      !! ** Action  : - computes the heat and freshwater/salt fluxes
91      !!              at the ice-ocean interface.
92      !!              - Update the ocean sbc
93      !!     
94      !! ** Outputs : - qsr     : sea heat flux:     solar
95      !!              - qns     : sea heat flux: non solar
96      !!              - emp     : freshwater budget: volume flux
97      !!              - sfx     : salt flux
98      !!              - fr_i    : ice fraction
99      !!              - tn_ice  : sea-ice surface temperature
100      !!              - alb_ice : sea-ice albedo (recomputed only for coupled mode)
101      !!
102      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
103      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
104      !!              These refs are now obsolete since everything has been revised
105      !!              The ref should be Rousset et al., 2015
106      !!---------------------------------------------------------------------
107      INTEGER, INTENT(in) ::   kt   ! number of iteration
108      !
109      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
110      REAL(wp) ::   zqmass           ! Heat flux associated with mass exchange ice->ocean (W.m-2)
111      REAL(wp) ::   zqsr             ! New solar flux received by the ocean
112      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_cs, zalb_os     ! 3D workspace
113      !!---------------------------------------------------------------------
114      IF( nn_timing == 1 )  CALL timing_start('ice_update_flx')
115
116      IF( kt == nit000 .AND. lwp ) THEN
117         WRITE(numout,*)
118         WRITE(numout,*)'ice_update_flx'
119         WRITE(numout,*)'~~~~~~~~~~~~~~'
120      ENDIF
121
122      ! --- case we bypass ice thermodynamics --- !
123      IF( .NOT. ln_limthd ) THEN   ! we suppose ice is impermeable => ocean is isolated from atmosphere
124         hfx_in   (:,:)   = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
125         hfx_out  (:,:)   = ( 1._wp - at_i_b(:,:) ) *   qns_oce(:,:)                  + qemp_oce(:,:)
126         ftr_ice  (:,:,:) = 0._wp
127         emp_ice  (:,:)   = 0._wp
128         qemp_ice (:,:)   = 0._wp
129         qevap_ice(:,:,:) = 0._wp
130      ENDIF
131     
132      DO jj = 1, jpj
133         DO ji = 1, jpi
134
135            !------------------------------------------!
136            !      heat flux at the ocean surface      !
137            !------------------------------------------!
138            ! Solar heat flux reaching the ocean = zqsr (W.m-2)
139            !---------------------------------------------------
140            zqsr = qsr_tot(ji,jj)
141            DO jl = 1, jpl
142               zqsr = zqsr - a_i_b(ji,jj,jl) * (  qsr_ice(ji,jj,jl) - ftr_ice(ji,jj,jl) ) 
143            END DO
144!!gm  why not like almost everywhere else :
145!!gm        zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * (  qsr_ice(ji,jj,:) - ftr_ice(ji,jj,:) )
146
147            ! Total heat flux reaching the ocean = hfx_out (W.m-2)
148            !---------------------------------------------------
149            zqmass         = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC)
150            hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr
151
152            ! Add the residual from heat diffusion equation and sublimation (W.m-2)
153            !----------------------------------------------------------------------
154            hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj) +   &
155               &           ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) )
156
157            ! New qsr and qns used to compute the oceanic heat flux at the next time step
158            !----------------------------------------------------------------------------
159            qsr(ji,jj) = zqsr                                     
160            qns(ji,jj) = hfx_out(ji,jj) - zqsr             
161
162            ! Mass flux at the atm. surface       
163            !-----------------------------------
164            wfx_sub(ji,jj) = wfx_snw_sub(ji,jj) + wfx_ice_sub(ji,jj)
165
166            ! Mass flux at the ocean surface     
167            !------------------------------------
168            !  case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED)
169            !  -------------------------------------------------------------------------------------
170            !  The idea of this approach is that the system that we consider is the ICE-OCEAN system
171            !  Thus  FW  flux  =  External ( E-P+snow melt)
172            !       Salt flux  =  Exchanges in the ice-ocean system then converted into FW
173            !                     Associated to Ice formation AND Ice melting
174            !                     Even if i see Ice melting as a FW and SALT flux
175            !       
176            ! mass flux from ice/ocean
177            wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   &
178               &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) 
179
180            IF ( ln_pnd_fw )   wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj)
181
182            ! add the snow melt water to snow mass flux to the ocean
183            wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj)
184
185            ! mass flux at the ocean/ice interface
186            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
187            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)
188
189
190            ! Salt flux at the ocean surface     
191            !------------------------------------------
192            sfx(ji,jj) = sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj)   &
193               &       + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj)
194           
195            ! Mass of snow and ice per unit area   
196            !----------------------------------------
197            snwice_mass_b(ji,jj) = snwice_mass(ji,jj)       ! save mass from the previous ice time step
198            !                                               ! new mass per unit area
199            snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj)  ) 
200            !                                               ! time evolution of snow+ice mass
201            snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_rdtice
202           
203         END DO
204      END DO
205
206      !-----------------------------------------------!
207      !   Storing the transmitted variables           !
208      !-----------------------------------------------!
209      fr_i  (:,:)   = at_i(:,:)             ! Sea-ice fraction           
210      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                     
211
212      !------------------------------------------------------------------------!
213      !    Snow/ice albedo (only if sent to coupler, useless in forced mode)   !
214      !------------------------------------------------------------------------!
215      CALL ice_alb( t_su, ht_i, ht_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
216      !
217      alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
218
219      !                    ! conservation test
220      IF( ln_limdiachk .AND. .NOT. ln_bdy)   CALL ice_cons_final( 'iceupdate' )
221      !                    ! control prints
222      IF( ln_limctl )   CALL ice_prt( kt, iiceprt, jiceprt, 3, ' - Final state ice_update - ' )
223      IF( ln_ctl    )   CALL ice_prt3D( 'iceupdate' )
224      !
225      IF( nn_timing == 1 )   CALL timing_stop('ice_update_flx')
226      !
227   END SUBROUTINE ice_update_flx
228
229
230   SUBROUTINE ice_update_tau( kt, pu_oce, pv_oce )
231      !!-------------------------------------------------------------------
232      !!                ***  ROUTINE ice_update_tau ***
233      !! 
234      !! ** Purpose : Update the ocean surface stresses due to the ice
235      !!         
236      !! ** Action  : * at each ice time step (every nn_fsbc time step):
237      !!                - compute the modulus of ice-ocean relative velocity
238      !!                  (*rho*Cd) at T-point (C-grid) or I-point (B-grid)
239      !!                      tmod_io = rhoco * | U_ice-U_oce |
240      !!                - update the modulus of stress at ocean surface
241      !!                      taum = (1-a) * taum + a * tmod_io * | U_ice-U_oce |
242      !!              * at each ocean time step (every kt):
243      !!                  compute linearized ice-ocean stresses as
244      !!                      Utau = tmod_io * | U_ice - pU_oce |
245      !!                using instantaneous current ocean velocity (usually before)
246      !!
247      !!    NB: - ice-ocean rotation angle no more allowed
248      !!        - here we make an approximation: taum is only computed every ice time step
249      !!          This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
250      !!          to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
251      !!
252      !! ** Outputs : - utau, vtau   : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
253      !!              - taum         : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
254      !!---------------------------------------------------------------------
255      INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index
256      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents
257      !
258      INTEGER  ::   ji, jj   ! dummy loop indices
259      REAL(wp) ::   zat_u, zutau_ice, zu_t, zmodt   ! local scalar
260      REAL(wp) ::   zat_v, zvtau_ice, zv_t, zrhoco  !   -      -
261      !!---------------------------------------------------------------------
262
263      IF( nn_timing == 1 )  CALL timing_start('ice_update_tau')
264
265      IF( kt == nit000 .AND. lwp ) THEN
266         WRITE(numout,*)
267         WRITE(numout,*)'ice_update_tau'
268         WRITE(numout,*)'~~~~~~~~~~~~~~'
269      ENDIF
270
271      zrhoco = rau0 * rn_cio
272      !
273      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
274         DO jj = 2, jpjm1                             !* update the modulus of stress at ocean surface (T-point)
275            DO ji = fs_2, fs_jpim1
276               !                                               ! 2*(U_ice-U_oce) at T-point
277               zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)   
278               zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) 
279               !                                              ! |U_ice-U_oce|^2
280               zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  )
281               !                                               ! update the ocean stress modulus
282               taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * zrhoco * zmodt
283               tmod_io(ji,jj) = zrhoco * SQRT( zmodt )          ! rhoco * |U_ice-U_oce| at T-point
284            END DO
285         END DO
286         CALL lbc_lnk_multi( taum, 'T', 1., tmod_io, 'T', 1. )
287         !
288         utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step
289         vtau_oce(:,:) = vtau(:,:)
290         !
291      ENDIF
292      !
293      !                                      !==  every ocean time-step  ==!
294      !
295      DO jj = 2, jpjm1                                !* update the stress WITHOUT a ice-ocean rotation angle
296         DO ji = fs_2, fs_jpim1   ! Vect. Opt.
297            zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp   ! ice area at u and V-points
298            zat_v  = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5_wp
299            !                                                   ! linearized quadratic drag formulation
300            zutau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
301            zvtau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
302            !                                                   ! stresses at the ocean surface
303            utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
304            vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
305         END DO
306      END DO
307      CALL lbc_lnk_multi( utau, 'U', -1., vtau, 'V', -1. )   ! lateral boundary condition
308      !
309      IF( nn_timing == 1 )  CALL timing_stop('ice_update_tau')
310     
311   END SUBROUTINE ice_update_tau
312
313
314   SUBROUTINE ice_update_init
315      !!-------------------------------------------------------------------
316      !!                  ***  ROUTINE ice_update_init  ***
317      !!             
318      !! ** Purpose :   ???
319      !!
320      !!-------------------------------------------------------------------
321      INTEGER  ::   ji, jj, jk               ! dummy loop indices
322      REAL(wp) ::   zcoefu, zcoefv, zcoeff   ! local scalar
323      !!-------------------------------------------------------------------
324      !
325      IF(lwp) WRITE(numout,*)
326      IF(lwp) WRITE(numout,*) 'ice_update_init :   sea-ice   ????'
327      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~   '
328
329      !                                      ! allocate ice_update array
330      IF( ice_update_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_update_init : unable to allocate standard arrays' )
331      !
332      soce_0(:,:) = soce                     ! constant SSS and ice salinity used in levitating case 0 (i.e. virtual salt flux)
333      sice_0(:,:) = sice
334      WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   &   ! reduced values in the Baltic Sea area
335         &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         ) 
336         soce_0(:,:) = 4._wp
337         sice_0(:,:) = 2._wp
338      END WHERE
339      !
340      IF( .NOT.ln_rstart ) THEN              ! set
341         !
342         snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )   ! snow+ice mass
343         snwice_mass_b(:,:) = snwice_mass(:,:)
344         !
345         IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
346            sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
347            sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
348
349!!gm I really don't like this stuff here...  Find a way to put that elsewhere or differently
350!!gm
351            IF( .NOT.ln_linssh ) THEN
352               DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
353                  e3t_n(:,:,jk) = e3t_0(:,:,jk) * (  1._wp + sshn(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1) )  )
354                  e3t_b(:,:,jk) = e3t_0(:,:,jk) * (  1._wp + sshb(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1) )  )
355               END DO
356               e3t_a(:,:,:) = e3t_b(:,:,:)
357!!gm  we are in no-restart case, so sshn=sshb   ==>> faster calculation:
358!!    REAL(wp) ::   ze3t   ! local scalar
359!!    REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! workspace
360!!
361!!             WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
362!!             ELSEWHERE                ;   z2d(:,:) = 1._wp
363!!             END WHERE
364!!             DO jk = 1,jpkm1                     ! adjust initial vertical scale factors               
365!!                ze3t = e3t_0(:,:,jk) * z2d(:,:)
366!!                e3t_n(:,:,jk) = ze3t
367!!                e3t_b(:,:,jk) = ze3t
368!!                e3t_a(:,:,jk) = ze3t
369!!             END DO
370!!gm  but since it is only done at the initialisation....  just the following can be acceptable:
371!               DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
372!                  e3t_n(:,:,jk) = e3t_0(:,:,jk) * (  1._wp + sshn(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1))  )
373!               END DO
374!               e3t_b(:,:,:) = e3t_n(:,:,:)
375!               e3t_a(:,:,:) = e3t_n(:,:,:)
376!!gm end               
377               ! Reconstruction of all vertical scale factors at now and before time-steps
378               ! =========================================================================
379               ! Horizontal scale factor interpolations
380               ! --------------------------------------
381               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
382               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
383               CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
384               CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
385               CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
386               ! Vertical scale factor interpolations
387                 ! ------------------------------------
388               CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
389               CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
390               CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
391               CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
392               CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
393               ! t- and w- points depth
394               ! ----------------------
395!!gm not sure of that....
396               gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
397               gdepw_n(:,:,1) = 0.0_wp
398               gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
399               DO jk = 2, jpk
400                  gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk  )
401                  gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
402                  gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn (:,:)
403               END DO
404            ENDIF
405         ENDIF
406      ENDIF ! .NOT. ln_rstart
407      !
408   END SUBROUTINE ice_update_init
409
410#else
411   !!----------------------------------------------------------------------
412   !!   Default option         Dummy module          NO  LIM3 sea-ice model
413   !!----------------------------------------------------------------------
414#endif 
415
416   !!======================================================================
417END MODULE iceupdate
Note: See TracBrowser for help on using the repository browser.