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

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90 @ 9015

Last change on this file since 9015 was 9015, checked in by acc, 6 years ago

Branch dev_CNRS_2017. Tidy up timing calls. Correct a few inconsistent labels and fix incorrect call in trc_sub_stp. Remove timing from icethd_pnd which is not always called by all processors

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