source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE/iceupdate.F90 @ 13140

Last change on this file since 13140 was 12546, checked in by orioltp, 12 months ago

Adding precision specification in hardcoded reals and other modifications to allow compilation without forcing reals without precision specification to a certain value through compiler flags

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