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

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90 @ 9119

Last change on this file since 9119 was 9071, checked in by clem, 6 years ago

mask some outputs

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