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

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

Last change on this file since 8512 was 8512, checked in by clem, 7 years ago

changes in style - part5 - reaching the end

File size: 22.1 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: update fluxes (mass, salt and heat) at the ice-ocean interface'
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            ! Solar heat flux reaching the ocean = zqsr (W.m-2)
136            !---------------------------------------------------
137            zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - ftr_ice(ji,jj,:) ) )
138
139            ! Total heat flux reaching the ocean = hfx_out (W.m-2)
140            !---------------------------------------------------
141            zqmass         = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC)
142            hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr
143
144            ! Add the residual from heat diffusion equation and sublimation (W.m-2)
145            !----------------------------------------------------------------------
146            hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj) +   &
147               &           ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) )
148
149            ! New qsr and qns used to compute the oceanic heat flux at the next time step
150            !----------------------------------------------------------------------------
151            qsr(ji,jj) = zqsr                                     
152            qns(ji,jj) = hfx_out(ji,jj) - zqsr             
153
154            ! Mass flux at the atm. surface       
155            !-----------------------------------
156            wfx_sub(ji,jj) = wfx_snw_sub(ji,jj) + wfx_ice_sub(ji,jj)
157
158            ! Mass flux at the ocean surface     
159            !------------------------------------
160            !  case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED)
161            !  -------------------------------------------------------------------------------------
162            !  The idea of this approach is that the system that we consider is the ICE-OCEAN system
163            !  Thus  FW  flux  =  External ( E-P+snow melt)
164            !       Salt flux  =  Exchanges in the ice-ocean system then converted into FW
165            !                     Associated to Ice formation AND Ice melting
166            !                     Even if i see Ice melting as a FW and SALT flux
167            !       
168            ! mass flux from ice/ocean
169            wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   &
170               &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) 
171
172            IF ( ln_pnd_fw )   wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj)
173
174            ! add the snow melt water to snow mass flux to the ocean
175            wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj)
176
177            ! mass flux at the ocean/ice interface
178            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
179            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)
180
181
182            ! Salt flux at the ocean surface     
183            !------------------------------------------
184            sfx(ji,jj) = sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj)   &
185               &       + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj)
186           
187            ! Mass of snow and ice per unit area   
188            !----------------------------------------
189            snwice_mass_b(ji,jj) = snwice_mass(ji,jj)       ! save mass from the previous ice time step
190            !                                               ! new mass per unit area
191            snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj)  ) 
192            !                                               ! time evolution of snow+ice mass
193            snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_rdtice
194           
195         END DO
196      END DO
197
198      ! Storing the transmitted variables
199      !----------------------------------
200      fr_i  (:,:)   = at_i(:,:)             ! Sea-ice fraction           
201      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                     
202
203      ! Snow/ice albedo (only if sent to coupler, useless in forced mode)
204      !------------------------------------------------------------------
205      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
206      !
207      alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
208
209      !                    ! conservation test
210      IF( ln_limdiachk .AND. .NOT. ln_bdy)   CALL ice_cons_final( 'iceupdate' )
211      !                    ! control prints
212      IF( ln_limctl )   CALL ice_prt( kt, iiceprt, jiceprt, 3, ' - Final state ice_update - ' )
213      IF( ln_ctl    )   CALL ice_prt3D( 'iceupdate' )
214      !
215      IF( nn_timing == 1 )   CALL timing_stop('ice_update_flx')
216      !
217   END SUBROUTINE ice_update_flx
218
219
220   SUBROUTINE ice_update_tau( kt, pu_oce, pv_oce )
221      !!-------------------------------------------------------------------
222      !!                ***  ROUTINE ice_update_tau ***
223      !! 
224      !! ** Purpose : Update the ocean surface stresses due to the ice
225      !!         
226      !! ** Action  : * at each ice time step (every nn_fsbc time step):
227      !!                - compute the modulus of ice-ocean relative velocity
228      !!                  (*rho*Cd) at T-point (C-grid) or I-point (B-grid)
229      !!                      tmod_io = rhoco * | U_ice-U_oce |
230      !!                - update the modulus of stress at ocean surface
231      !!                      taum = (1-a) * taum + a * tmod_io * | U_ice-U_oce |
232      !!              * at each ocean time step (every kt):
233      !!                  compute linearized ice-ocean stresses as
234      !!                      Utau = tmod_io * | U_ice - pU_oce |
235      !!                using instantaneous current ocean velocity (usually before)
236      !!
237      !!    NB: - ice-ocean rotation angle no more allowed
238      !!        - here we make an approximation: taum is only computed every ice time step
239      !!          This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
240      !!          to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
241      !!
242      !! ** Outputs : - utau, vtau   : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
243      !!              - taum         : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
244      !!---------------------------------------------------------------------
245      INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index
246      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents
247      !
248      INTEGER  ::   ji, jj   ! dummy loop indices
249      REAL(wp) ::   zat_u, zutau_ice, zu_t, zmodt   ! local scalar
250      REAL(wp) ::   zat_v, zvtau_ice, zv_t, zrhoco  !   -      -
251      !!---------------------------------------------------------------------
252
253      IF( nn_timing == 1 )  CALL timing_start('ice_update_tau')
254
255      IF( kt == nit000 .AND. lwp ) THEN
256         WRITE(numout,*)
257         WRITE(numout,*)'ice_update_tau: update stress at the ice-ocean interface'
258         WRITE(numout,*)'~~~~~~~~~~~~~~'
259      ENDIF
260
261      zrhoco = rau0 * rn_cio
262      !
263      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
264         DO jj = 2, jpjm1                             !* update the modulus of stress at ocean surface (T-point)
265            DO ji = fs_2, fs_jpim1
266               !                                               ! 2*(U_ice-U_oce) at T-point
267               zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)   
268               zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) 
269               !                                              ! |U_ice-U_oce|^2
270               zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  )
271               !                                               ! update the ocean stress modulus
272               taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * zrhoco * zmodt
273               tmod_io(ji,jj) = zrhoco * SQRT( zmodt )          ! rhoco * |U_ice-U_oce| at T-point
274            END DO
275         END DO
276         CALL lbc_lnk_multi( taum, 'T', 1., tmod_io, 'T', 1. )
277         !
278         utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step
279         vtau_oce(:,:) = vtau(:,:)
280         !
281      ENDIF
282      !
283      !                                      !==  every ocean time-step  ==!
284      !
285      DO jj = 2, jpjm1                                !* update the stress WITHOUT a ice-ocean rotation angle
286         DO ji = fs_2, fs_jpim1   ! Vect. Opt.
287            zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp   ! ice area at u and V-points
288            zat_v  = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5_wp
289            !                                                   ! linearized quadratic drag formulation
290            zutau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
291            zvtau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
292            !                                                   ! stresses at the ocean surface
293            utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
294            vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
295         END DO
296      END DO
297      CALL lbc_lnk_multi( utau, 'U', -1., vtau, 'V', -1. )   ! lateral boundary condition
298      !
299      IF( nn_timing == 1 )  CALL timing_stop('ice_update_tau')
300     
301   END SUBROUTINE ice_update_tau
302
303
304   SUBROUTINE ice_update_init
305      !!-------------------------------------------------------------------
306      !!                  ***  ROUTINE ice_update_init  ***
307      !!             
308      !! ** Purpose :   ???
309      !!
310      !!-------------------------------------------------------------------
311      INTEGER  ::   ji, jj, jk               ! dummy loop indices
312      REAL(wp) ::   zcoefu, zcoefv, zcoeff   ! local scalar
313      !!-------------------------------------------------------------------
314      !
315      IF(lwp) WRITE(numout,*)
316      IF(lwp) WRITE(numout,*) 'ice_update_init: ???? '
317      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
318
319      !                                      ! allocate ice_update array
320      IF( ice_update_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_update_init : unable to allocate standard arrays' )
321      !
322      soce_0(:,:) = soce                     ! constant SSS and ice salinity used in levitating case 0 (i.e. virtual salt flux)
323      sice_0(:,:) = sice
324      WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   &   ! reduced values in the Baltic Sea area
325         &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         ) 
326         soce_0(:,:) = 4._wp
327         sice_0(:,:) = 2._wp
328      END WHERE
329      !
330      IF( .NOT.ln_rstart ) THEN              ! set
331         !
332         snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )   ! snow+ice mass
333         snwice_mass_b(:,:) = snwice_mass(:,:)
334         !
335         IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
336            sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
337            sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
338
339!!gm I really don't like this stuff here...  Find a way to put that elsewhere or differently
340!!gm
341            IF( .NOT.ln_linssh ) THEN
342               DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
343                  e3t_n(:,:,jk) = e3t_0(:,:,jk) * (  1._wp + sshn(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1) )  )
344                  e3t_b(:,:,jk) = e3t_0(:,:,jk) * (  1._wp + sshb(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1) )  )
345               END DO
346               e3t_a(:,:,:) = e3t_b(:,:,:)
347!!gm  we are in no-restart case, so sshn=sshb   ==>> faster calculation:
348!!    REAL(wp) ::   ze3t   ! local scalar
349!!    REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! workspace
350!!
351!!             WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
352!!             ELSEWHERE                ;   z2d(:,:) = 1._wp
353!!             END WHERE
354!!             DO jk = 1,jpkm1                     ! adjust initial vertical scale factors               
355!!                ze3t = e3t_0(:,:,jk) * z2d(:,:)
356!!                e3t_n(:,:,jk) = ze3t
357!!                e3t_b(:,:,jk) = ze3t
358!!                e3t_a(:,:,jk) = ze3t
359!!             END DO
360!!gm  but since it is only done at the initialisation....  just the following can be acceptable:
361!               DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
362!                  e3t_n(:,:,jk) = e3t_0(:,:,jk) * (  1._wp + sshn(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1))  )
363!               END DO
364!               e3t_b(:,:,:) = e3t_n(:,:,:)
365!               e3t_a(:,:,:) = e3t_n(:,:,:)
366!!gm end               
367               ! Reconstruction of all vertical scale factors at now and before time-steps
368               ! =========================================================================
369               ! Horizontal scale factor interpolations
370               ! --------------------------------------
371               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
372               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
373               CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
374               CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
375               CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
376               ! Vertical scale factor interpolations
377                 ! ------------------------------------
378               CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
379               CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
380               CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
381               CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
382               CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
383               ! t- and w- points depth
384               ! ----------------------
385!!gm not sure of that....
386               gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
387               gdepw_n(:,:,1) = 0.0_wp
388               gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
389               DO jk = 2, jpk
390                  gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk  )
391                  gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
392                  gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn (:,:)
393               END DO
394            ENDIF
395         ENDIF
396      ENDIF ! .NOT. ln_rstart
397      !
398   END SUBROUTINE ice_update_init
399
400#else
401   !!----------------------------------------------------------------------
402   !!   Default option         Dummy module          NO  LIM3 sea-ice model
403   !!----------------------------------------------------------------------
404#endif 
405
406   !!======================================================================
407END MODULE iceupdate
Note: See TracBrowser for help on using the repository browser.