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

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

changes in style - part6 - pure cosmetics

File size: 20.4 KB
Line 
1MODULE iceupdate
2   !!======================================================================
3   !!                       ***  MODULE iceupdate   ***
4   !!  Sea-ice :   computation of the flux at the sea ice/ocean interface
5   !!======================================================================
6   !! History :   -   ! 2006-07 (M. Vancoppelle)  LIM3 original code
7   !!            3.0  ! 2008-03 (C. Tallandier)  surface module
8   !!             -   ! 2008-04 (C. Tallandier)  split in 2 + new ice-ocean coupling
9   !!            3.3  ! 2010-05 (G. Madec) decrease ocean & ice reference salinities in the Baltic sea
10   !!                 !                  + simplification of the ice-ocean stress calculation
11   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
12   !!             -   ! 2012    (D. Iovino) salt flux change
13   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
14   !!            3.5  ! 2012-10 (A. Coward, G. Madec) salt fluxes ; ice+snow mass
15   !!----------------------------------------------------------------------
16#if defined key_lim3
17   !!----------------------------------------------------------------------
18   !!   'key_lim3'                                       ESIM sea-ice model
19   !!----------------------------------------------------------------------
20   !!   ice_update_alloc : allocate the iceupdate arrays
21   !!   ice_update_init  : initialisation
22   !!   ice_update_flx   : updates mass, heat and salt fluxes at the ocean surface
23   !!   ice_update_tau   : update i- and j-stresses, and its modulus at the ocean surface
24   !!----------------------------------------------------------------------
25   USE oce     , ONLY : sshn, sshb
26   USE phycst         ! physical constants
27   USE dom_oce        ! ocean domain
28   USE ice            ! sea-ice: variables
29!!gm  It should be probably better to pass these variable in argument of the routine,
30!!gm  rather than having this long list in USE. This will also highlight what is updated, and what is just use.
31   USE sbc_ice , ONLY : emp_oce, qns_oce, qsr_oce , qemp_oce ,                             &
32      &                 emp_ice, qsr_ice, qemp_ice, qevap_ice, alb_ice, tn_ice, cldf_ice,  &
33      &                 snwice_mass, snwice_mass_b, snwice_fmass
34   USE sbc_oce , ONLY : nn_fsbc, ln_ice_embd, sfx, fr_i, qsr_tot, qns, qsr, fmmflx, emp, taum, utau, vtau
35!!gm end
36   USE sbccpl         ! Surface boundary condition: coupled interface
37   USE icealb         ! sea-ice: albedo parameters
38   USE traqsr         ! add penetration of solar flux in the calculation of heat budget
39   USE icectl         ! sea-ice: control prints
40   USE bdy_oce , ONLY : ln_bdy
41   !
42   USE in_out_manager ! I/O manager
43   USE iom            ! I/O manager library
44   USE lib_mpp        ! MPP library
45   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
46   USE lbclnk         ! lateral boundary conditions (or mpp links)
47   USE timing         ! Timing
48
49   IMPLICIT NONE
50   PRIVATE
51
52   PUBLIC   ice_update_init   ! called by ice_init
53   PUBLIC   ice_update_flx    ! called by ice_stp
54   PUBLIC   ice_update_tau    ! called by ice_stp
55
56   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress     [N/m2]
57   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmod_io              ! modulus of the ice-ocean velocity   [m/s]
58
59   !! * Substitutions
60#  include "vectopt_loop_substitute.h90"
61   !!----------------------------------------------------------------------
62   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
63   !! $Id: iceupdate.F90 8411 2017-08-07 16:09:12Z clem $
64   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
65   !!----------------------------------------------------------------------
66CONTAINS
67
68   INTEGER FUNCTION ice_update_alloc()
69      !!-------------------------------------------------------------------
70      !!             ***  ROUTINE ice_update_alloc ***
71      !!-------------------------------------------------------------------
72      ALLOCATE( utau_oce(jpi,jpj), vtau_oce(jpi,jpj), tmod_io(jpi,jpj), STAT=ice_update_alloc )
73      !
74      IF( lk_mpp                )   CALL mpp_sum( ice_update_alloc )
75      IF( ice_update_alloc /= 0 )   CALL ctl_warn('ice_update_alloc: failed to allocate arrays')
76      !
77   END FUNCTION ice_update_alloc
78
79
80   SUBROUTINE ice_update_flx( kt )
81      !!-------------------------------------------------------------------
82      !!                ***  ROUTINE ice_update_flx ***
83      !! 
84      !! ** Purpose :   Update the surface ocean boundary condition for heat
85      !!                salt and mass over areas where sea-ice is non-zero
86      !!         
87      !! ** Action  : - computes the heat and freshwater/salt fluxes
88      !!                at the ice-ocean interface.
89      !!              - Update the ocean sbc
90      !!     
91      !! ** Outputs : - qsr     : sea heat flux:     solar
92      !!              - qns     : sea heat flux: non solar
93      !!              - emp     : freshwater budget: volume flux
94      !!              - sfx     : salt flux
95      !!              - fr_i    : ice fraction
96      !!              - tn_ice  : sea-ice surface temperature
97      !!              - alb_ice : sea-ice albedo (recomputed only for coupled mode)
98      !!
99      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
100      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
101      !!              These refs are now obsolete since everything has been revised
102      !!              The ref should be Rousset et al., 2015
103      !!---------------------------------------------------------------------
104      INTEGER, INTENT(in) ::   kt   ! number of iteration
105      !
106      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
107      REAL(wp) ::   zqmass           ! Heat flux associated with mass exchange ice->ocean (W.m-2)
108      REAL(wp) ::   zqsr             ! New solar flux received by the ocean
109      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_cs, zalb_os     ! 3D workspace
110      !!---------------------------------------------------------------------
111      IF( nn_timing == 1 )  CALL timing_start('ice_update')
112
113      IF( kt == nit000 .AND. lwp ) THEN
114         WRITE(numout,*)
115         WRITE(numout,*)'ice_update_flx: update fluxes (mass, salt and heat) at the ice-ocean interface'
116         WRITE(numout,*)'~~~~~~~~~~~~~~'
117      ENDIF
118
119      ! --- case we bypass ice thermodynamics --- !
120      IF( .NOT. ln_icethd ) THEN   ! we suppose ice is impermeable => ocean is isolated from atmosphere
121         hfx_in   (:,:)   = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
122         hfx_out  (:,:)   = ( 1._wp - at_i_b(:,:) ) *   qns_oce(:,:)                  + qemp_oce(:,:)
123         ftr_ice  (:,:,:) = 0._wp
124         emp_ice  (:,:)   = 0._wp
125         qemp_ice (:,:)   = 0._wp
126         qevap_ice(:,:,:) = 0._wp
127      ENDIF
128     
129      DO jj = 1, jpj
130         DO ji = 1, jpi
131
132            ! Solar heat flux reaching the ocean = zqsr (W.m-2)
133            !---------------------------------------------------
134            zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - ftr_ice(ji,jj,:) ) )
135
136            ! Total heat flux reaching the ocean = hfx_out (W.m-2)
137            !---------------------------------------------------
138            zqmass         = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC)
139            hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr
140
141            ! Add the residual from heat diffusion equation and sublimation (W.m-2)
142            !----------------------------------------------------------------------
143            hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj) +   &
144               &           ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) )
145
146            ! New qsr and qns used to compute the oceanic heat flux at the next time step
147            !----------------------------------------------------------------------------
148            qsr(ji,jj) = zqsr                                     
149            qns(ji,jj) = hfx_out(ji,jj) - zqsr             
150
151            ! Mass flux at the atm. surface       
152            !-----------------------------------
153            wfx_sub(ji,jj) = wfx_snw_sub(ji,jj) + wfx_ice_sub(ji,jj)
154
155            ! Mass flux at the ocean surface     
156            !------------------------------------
157            !  case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED)
158            !  -------------------------------------------------------------------------------------
159            !  The idea of this approach is that the system that we consider is the ICE-OCEAN system
160            !  Thus  FW  flux  =  External ( E-P+snow melt)
161            !       Salt flux  =  Exchanges in the ice-ocean system then converted into FW
162            !                     Associated to Ice formation AND Ice melting
163            !                     Even if i see Ice melting as a FW and SALT flux
164            !       
165            ! mass flux from ice/ocean
166            wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   &
167               &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) 
168
169            IF ( ln_pnd_fw )   wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj)
170
171            ! add the snow melt water to snow mass flux to the ocean
172            wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj)
173
174            ! mass flux at the ocean/ice interface
175            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
176            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)
177
178
179            ! Salt flux at the ocean surface     
180            !------------------------------------------
181            sfx(ji,jj) = sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj)   &
182               &       + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj)
183           
184            ! Mass of snow and ice per unit area   
185            !----------------------------------------
186            snwice_mass_b(ji,jj) = snwice_mass(ji,jj)       ! save mass from the previous ice time step
187            !                                               ! new mass per unit area
188            snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj)  ) 
189            !                                               ! time evolution of snow+ice mass
190            snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_rdtice
191           
192         END DO
193      END DO
194
195      ! Storing the transmitted variables
196      !----------------------------------
197      fr_i  (:,:)   = at_i(:,:)             ! Sea-ice fraction           
198      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                     
199
200      ! Snow/ice albedo (only if sent to coupler, useless in forced mode)
201      !------------------------------------------------------------------
202      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
203      !
204      alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
205      !
206      IF( lrst_ice ) THEN                       !* write snwice_mass fields in the restart file
207         CALL update_rst( 'WRITE', kt )
208      ENDIF
209      !
210      ! controls
211      IF( ln_icediachk .AND. .NOT. ln_bdy)   CALL ice_cons_final('iceupdate')                                       ! conservation
212      IF( ln_icectl                      )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints
213      IF( ln_ctl                         )   CALL ice_prt3D     ('iceupdate')                                       ! prints
214      IF( nn_timing == 1                 )   CALL timing_stop   ('ice_update')                                      ! timing
215      !
216   END SUBROUTINE ice_update_flx
217
218
219   SUBROUTINE ice_update_tau( kt, pu_oce, pv_oce )
220      !!-------------------------------------------------------------------
221      !!                ***  ROUTINE ice_update_tau ***
222      !! 
223      !! ** Purpose : Update the ocean surface stresses due to the ice
224      !!         
225      !! ** Action  : * at each ice time step (every nn_fsbc time step):
226      !!                - compute the modulus of ice-ocean relative velocity
227      !!                  (*rho*Cd) at T-point (C-grid) or I-point (B-grid)
228      !!                      tmod_io = rhoco * | U_ice-U_oce |
229      !!                - update the modulus of stress at ocean surface
230      !!                      taum = (1-a) * taum + a * tmod_io * | U_ice-U_oce |
231      !!              * at each ocean time step (every kt):
232      !!                  compute linearized ice-ocean stresses as
233      !!                      Utau = tmod_io * | U_ice - pU_oce |
234      !!                using instantaneous current ocean velocity (usually before)
235      !!
236      !!    NB: - ice-ocean rotation angle no more allowed
237      !!        - here we make an approximation: taum is only computed every ice time step
238      !!          This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
239      !!          to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
240      !!
241      !! ** Outputs : - utau, vtau   : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
242      !!              - taum         : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
243      !!---------------------------------------------------------------------
244      INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index
245      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents
246      !
247      INTEGER  ::   ji, jj   ! dummy loop indices
248      REAL(wp) ::   zat_u, zutau_ice, zu_t, zmodt   ! local scalar
249      REAL(wp) ::   zat_v, zvtau_ice, zv_t, zrhoco  !   -      -
250      !!---------------------------------------------------------------------
251
252      IF( nn_timing == 1 )  CALL timing_start('ice_update')
253
254      IF( kt == nit000 .AND. lwp ) THEN
255         WRITE(numout,*)
256         WRITE(numout,*)'ice_update_tau: update stress at the ice-ocean interface'
257         WRITE(numout,*)'~~~~~~~~~~~~~~'
258      ENDIF
259
260      zrhoco = rau0 * rn_cio
261      !
262      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
263         DO jj = 2, jpjm1                             !* update the modulus of stress at ocean surface (T-point)
264            DO ji = fs_2, fs_jpim1
265               !                                               ! 2*(U_ice-U_oce) at T-point
266               zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)   
267               zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) 
268               !                                              ! |U_ice-U_oce|^2
269               zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  )
270               !                                               ! update the ocean stress modulus
271               taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * zrhoco * zmodt
272               tmod_io(ji,jj) = zrhoco * SQRT( zmodt )          ! rhoco * |U_ice-U_oce| at T-point
273            END DO
274         END DO
275         CALL lbc_lnk_multi( taum, 'T', 1., tmod_io, 'T', 1. )
276         !
277         utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step
278         vtau_oce(:,:) = vtau(:,:)
279         !
280      ENDIF
281      !
282      !                                      !==  every ocean time-step  ==!
283      !
284      DO jj = 2, jpjm1                                !* update the stress WITHOUT a ice-ocean rotation angle
285         DO ji = fs_2, fs_jpim1   ! Vect. Opt.
286            zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp   ! ice area at u and V-points
287            zat_v  = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5_wp
288            !                                                   ! linearized quadratic drag formulation
289            zutau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
290            zvtau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
291            !                                                   ! stresses at the ocean surface
292            utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
293            vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
294         END DO
295      END DO
296      CALL lbc_lnk_multi( utau, 'U', -1., vtau, 'V', -1. )   ! lateral boundary condition
297      !
298      IF( nn_timing == 1 )  CALL timing_stop('ice_update')
299     
300   END SUBROUTINE ice_update_tau
301
302
303   SUBROUTINE ice_update_init
304      !!-------------------------------------------------------------------
305      !!                  ***  ROUTINE ice_update_init  ***
306      !!             
307      !! ** Purpose :   ???
308      !!
309      !!-------------------------------------------------------------------
310      INTEGER  ::   ji, jj, jk               ! dummy loop indices
311      REAL(wp) ::   zcoefu, zcoefv, zcoeff   ! local scalar
312      !!-------------------------------------------------------------------
313      !
314      IF(lwp) WRITE(numout,*)
315      IF(lwp) WRITE(numout,*) 'ice_update_init: ???? '
316      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
317
318      !                                      ! allocate ice_update array
319      IF( ice_update_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_update_init : unable to allocate standard arrays' )
320      !
321      CALL update_rst( 'READ' )  !* read or initialize all required files
322      !
323   END SUBROUTINE ice_update_init
324
325   SUBROUTINE update_rst( cdrw, kt )
326      !!---------------------------------------------------------------------
327      !!                   ***  ROUTINE rhg_evp_rst  ***
328      !!                     
329      !! ** Purpose :   Read or write RHG file in restart file
330      !!
331      !! ** Method  :   use of IOM library
332      !!----------------------------------------------------------------------
333      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
334      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
335      !
336      INTEGER  ::   iter   ! local integer
337      INTEGER  ::   id1    ! local integer
338      !!----------------------------------------------------------------------
339      !
340      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
341         !                                   ! ---------------
342         IF( ln_rstart ) THEN                   !* Read the restart file
343            !
344            id1 = iom_varid( numrir, 'snwice_mass' , ldstop = .FALSE. )
345            !
346            IF( id1 > 0 ) THEN                       ! fields exist
347               CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass'  , snwice_mass   )
348               CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b )
349            ELSE                                     ! start from rest
350               IF(lwp) WRITE(numout,*) '   ==>>   previous run without snow-ice mass output then set it'
351               snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
352               snwice_mass_b(:,:) = snwice_mass(:,:)
353            ENDIF
354         ELSE                                   !* Start from rest
355            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set the snow-ice mass'
356            snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
357            snwice_mass_b(:,:) = snwice_mass(:,:)
358         ENDIF
359         !
360      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
361         !                                   ! -------------------
362         IF(lwp) WRITE(numout,*) '---- update-rst ----'
363         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
364         !
365         CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass'  , snwice_mass   )
366         CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b )
367         !
368      ENDIF
369      !
370   END SUBROUTINE update_rst
371
372#else
373   !!----------------------------------------------------------------------
374   !!   Default option         Dummy module           NO ESIM sea-ice model
375   !!----------------------------------------------------------------------
376#endif 
377
378   !!======================================================================
379END MODULE iceupdate
Note: See TracBrowser for help on using the repository browser.