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 @ 8531

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

changes in style - part6 - commits of the day

File size: 20.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         ! sea-ice: albedo parameters
39   USE traqsr         ! add penetration of solar flux in the calculation of heat budget
40   USE icectl         ! sea-ice: control prints
41   USE bdy_oce , ONLY : ln_bdy
42   !
43   USE in_out_manager ! I/O manager
44   USE iom            ! xIO server
45   USE lbclnk         ! ocean lateral boundary condition - MPP exchanges
46   USE lib_mpp        ! MPP library
47   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
48   USE timing         ! Timing
49
50   IMPLICIT NONE
51   PRIVATE
52
53   PUBLIC   ice_update_init   ! called by ice_init
54   PUBLIC   ice_update_flx    ! called by ice_stp
55   PUBLIC   ice_update_tau    ! called by ice_stp
56
57   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress     [N/m2]
58   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmod_io              ! modulus of the ice-ocean velocity   [m/s]
59
60   !! * Substitutions
61#  include "vectopt_loop_substitute.h90"
62   !!----------------------------------------------------------------------
63   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
64   !! $Id: iceupdate.F90 8411 2017-08-07 16:09:12Z clem $
65   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
66   !!----------------------------------------------------------------------
67CONTAINS
68
69   INTEGER FUNCTION ice_update_alloc()
70      !!-------------------------------------------------------------------
71      !!             ***  ROUTINE ice_update_alloc ***
72      !!-------------------------------------------------------------------
73      ALLOCATE( utau_oce(jpi,jpj), vtau_oce(jpi,jpj), tmod_io(jpi,jpj), STAT=ice_update_alloc )
74      !
75      IF( lk_mpp                )   CALL mpp_sum( ice_update_alloc )
76      IF( ice_update_alloc /= 0 )   CALL ctl_warn('ice_update_alloc: failed to allocate arrays')
77      !
78   END FUNCTION ice_update_alloc
79
80
81   SUBROUTINE ice_update_flx( kt )
82      !!-------------------------------------------------------------------
83      !!                ***  ROUTINE ice_update_flx ***
84      !! 
85      !! ** Purpose :   Update the surface ocean boundary condition for heat
86      !!                salt and mass over areas where sea-ice is non-zero
87      !!         
88      !! ** Action  : - computes the heat and freshwater/salt fluxes
89      !!                at the ice-ocean interface.
90      !!              - Update the ocean sbc
91      !!     
92      !! ** Outputs : - qsr     : sea heat flux:     solar
93      !!              - qns     : sea heat flux: non solar
94      !!              - emp     : freshwater budget: volume flux
95      !!              - sfx     : salt flux
96      !!              - fr_i    : ice fraction
97      !!              - tn_ice  : sea-ice surface temperature
98      !!              - alb_ice : sea-ice albedo (recomputed only for coupled mode)
99      !!
100      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
101      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
102      !!              These refs are now obsolete since everything has been revised
103      !!              The ref should be Rousset et al., 2015
104      !!---------------------------------------------------------------------
105      INTEGER, INTENT(in) ::   kt   ! number of iteration
106      !
107      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
108      REAL(wp) ::   zqmass           ! Heat flux associated with mass exchange ice->ocean (W.m-2)
109      REAL(wp) ::   zqsr             ! New solar flux received by the ocean
110      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_cs, zalb_os     ! 3D workspace
111      !!---------------------------------------------------------------------
112      IF( nn_timing == 1 )  CALL timing_start('ice_update_flx')
113
114      IF( kt == nit000 .AND. lwp ) THEN
115         WRITE(numout,*)
116         WRITE(numout,*)'ice_update_flx: update fluxes (mass, salt and heat) at the ice-ocean interface'
117         WRITE(numout,*)'~~~~~~~~~~~~~~'
118      ENDIF
119
120      ! --- case we bypass ice thermodynamics --- !
121      IF( .NOT. ln_icethd ) THEN   ! we suppose ice is impermeable => ocean is isolated from atmosphere
122         hfx_in   (:,:)   = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
123         hfx_out  (:,:)   = ( 1._wp - at_i_b(:,:) ) *   qns_oce(:,:)                  + qemp_oce(:,:)
124         ftr_ice  (:,:,:) = 0._wp
125         emp_ice  (:,:)   = 0._wp
126         qemp_ice (:,:)   = 0._wp
127         qevap_ice(:,:,:) = 0._wp
128      ENDIF
129     
130      DO jj = 1, jpj
131         DO ji = 1, jpi
132
133            ! Solar heat flux reaching the ocean = zqsr (W.m-2)
134            !---------------------------------------------------
135            zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - ftr_ice(ji,jj,:) ) )
136
137            ! Total heat flux reaching the ocean = hfx_out (W.m-2)
138            !---------------------------------------------------
139            zqmass         = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC)
140            hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr
141
142            ! Add the residual from heat diffusion equation and sublimation (W.m-2)
143            !----------------------------------------------------------------------
144            hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj) +   &
145               &           ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) )
146
147            ! New qsr and qns used to compute the oceanic heat flux at the next time step
148            !----------------------------------------------------------------------------
149            qsr(ji,jj) = zqsr                                     
150            qns(ji,jj) = hfx_out(ji,jj) - zqsr             
151
152            ! Mass flux at the atm. surface       
153            !-----------------------------------
154            wfx_sub(ji,jj) = wfx_snw_sub(ji,jj) + wfx_ice_sub(ji,jj)
155
156            ! Mass flux at the ocean surface     
157            !------------------------------------
158            !  case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED)
159            !  -------------------------------------------------------------------------------------
160            !  The idea of this approach is that the system that we consider is the ICE-OCEAN system
161            !  Thus  FW  flux  =  External ( E-P+snow melt)
162            !       Salt flux  =  Exchanges in the ice-ocean system then converted into FW
163            !                     Associated to Ice formation AND Ice melting
164            !                     Even if i see Ice melting as a FW and SALT flux
165            !       
166            ! mass flux from ice/ocean
167            wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   &
168               &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) 
169
170            IF ( ln_pnd_fw )   wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj)
171
172            ! add the snow melt water to snow mass flux to the ocean
173            wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj)
174
175            ! mass flux at the ocean/ice interface
176            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
177            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)
178
179
180            ! Salt flux at the ocean surface     
181            !------------------------------------------
182            sfx(ji,jj) = sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj)   &
183               &       + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj)
184           
185            ! Mass of snow and ice per unit area   
186            !----------------------------------------
187            snwice_mass_b(ji,jj) = snwice_mass(ji,jj)       ! save mass from the previous ice time step
188            !                                               ! new mass per unit area
189            snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj)  ) 
190            !                                               ! time evolution of snow+ice mass
191            snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_rdtice
192           
193         END DO
194      END DO
195
196      ! Storing the transmitted variables
197      !----------------------------------
198      fr_i  (:,:)   = at_i(:,:)             ! Sea-ice fraction           
199      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                     
200
201      ! Snow/ice albedo (only if sent to coupler, useless in forced mode)
202      !------------------------------------------------------------------
203      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
204      !
205      alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
206      !
207      IF( lrst_ice ) THEN                       !* write snwice_mass fields in the restart file
208         CALL update_rst( 'WRITE', kt )
209      ENDIF
210      !
211      ! controls
212      IF( ln_icediachk .AND. .NOT. ln_bdy)   CALL ice_cons_final('iceupdate')                                       ! conservation
213      IF( ln_icectl                      )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints
214      IF( ln_ctl                         )   CALL ice_prt3D     ('iceupdate')                                       ! prints
215      IF( nn_timing == 1                 )   CALL timing_stop   ('ice_update_flx')                                  ! timing
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      CALL update_rst( 'READ' )  !* read or initialize all required files
323      !
324   END SUBROUTINE ice_update_init
325
326   SUBROUTINE update_rst( cdrw, kt )
327      !!---------------------------------------------------------------------
328      !!                   ***  ROUTINE rhg_evp_rst  ***
329      !!                     
330      !! ** Purpose :   Read or write RHG file in restart file
331      !!
332      !! ** Method  :   use of IOM library
333      !!----------------------------------------------------------------------
334      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
335      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
336      !
337      INTEGER  ::   iter   ! local integer
338      INTEGER  ::   id1    ! local integer
339      !!----------------------------------------------------------------------
340      !
341      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
342         !                                   ! ---------------
343         IF( ln_rstart ) THEN                   !* Read the restart file
344            !
345            id1 = iom_varid( numrir, 'snwice_mass' , ldstop = .FALSE. )
346            !
347            IF( id1 > 0 ) THEN                       ! fields exist
348               CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass'  , snwice_mass   )
349               CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b )
350            ELSE                                     ! start from rest
351               IF(lwp) WRITE(numout,*) '   ==>>   previous run without snow-ice mass output then set it'
352               snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
353               snwice_mass_b(:,:) = snwice_mass(:,:)
354            ENDIF
355         ELSE                                   !* Start from rest
356            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set the snow-ice mass'
357            snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
358            snwice_mass_b(:,:) = snwice_mass(:,:)
359         ENDIF
360         !
361      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
362         !                                   ! -------------------
363         IF(lwp) WRITE(numout,*) '---- update-rst ----'
364         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
365         !
366         CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass'  , snwice_mass   )
367         CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b )
368         !
369      ENDIF
370      !
371   END SUBROUTINE update_rst
372
373#else
374   !!----------------------------------------------------------------------
375   !!   Default option         Dummy module          NO  LIM3 sea-ice model
376   !!----------------------------------------------------------------------
377#endif 
378
379   !!======================================================================
380END MODULE iceupdate
Note: See TracBrowser for help on using the repository browser.