Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
MODULE iceupdate
!!======================================================================
!! *** MODULE iceupdate ***
!! Sea-ice : computation of the flux at the sea ice/ocean interface
!!======================================================================
!! History : 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
!!----------------------------------------------------------------------
#if defined key_si3
!!----------------------------------------------------------------------
!! 'key_si3' SI3 sea-ice model
!!----------------------------------------------------------------------
!! ice_update_alloc : allocate the iceupdate arrays
!! ice_update_init : initialisation
!! ice_update_flx : updates mass, heat and salt fluxes at the ocean surface
!! ice_update_tau : update i- and j-stresses, and its modulus at the ocean surface
!!----------------------------------------------------------------------
USE phycst ! physical constants
USE dom_oce ! ocean domain
USE ice ! sea-ice: variables
USE sbc_ice ! Surface boundary condition: ice fields
USE sbc_oce ! Surface boundary condition: ocean fields
USE sbccpl ! Surface boundary condition: coupled interface
USE icealb ! sea-ice: albedo parameters
USE traqsr ! add penetration of solar flux in the calculation of heat budget
USE icectl ! sea-ice: control prints
USE zdfdrg , ONLY : ln_drgice_imp
!
USE in_out_manager ! I/O manager
USE iom ! I/O manager library
USE lib_mpp ! MPP library
USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
USE lbclnk ! lateral boundary conditions (or mpp links)
USE timing ! Timing
IMPLICIT NONE
PRIVATE
PUBLIC ice_update_init ! called by ice_init
PUBLIC ice_update_flx ! called by ice_stp
PUBLIC ice_update_tau ! called by ice_stp
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau_oce, vtau_oce ! air-ocean surface i- & j-stress [N/m2]
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmod_io ! modulus of the ice-ocean velocity [m/s]
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: iceupdate.F90 15385 2021-10-15 13:52:48Z clem $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
INTEGER FUNCTION ice_update_alloc()
!!-------------------------------------------------------------------
!! *** ROUTINE ice_update_alloc ***
!!-------------------------------------------------------------------
ALLOCATE( utau_oce(A2D(0)), vtau_oce(A2D(0)), tmod_io(jpi,jpj), STAT=ice_update_alloc )
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
!
CALL mpp_sum( 'iceupdate', ice_update_alloc )
IF( ice_update_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_update_alloc: failed to allocate arrays' )
!
END FUNCTION ice_update_alloc
SUBROUTINE ice_update_flx( kt )
!!-------------------------------------------------------------------
!! *** ROUTINE ice_update_flx ***
!!
!! ** Purpose : Update the surface ocean boundary condition for heat
!! salt and mass over areas where sea-ice is non-zero
!!
!! ** Action : - computes the heat and freshwater/salt fluxes
!! at the ice-ocean interface.
!! - Update the ocean sbc
!!
!! ** Outputs : - qsr : sea heat flux: solar
!! - qns : sea heat flux: non solar
!! - emp : freshwater budget: volume flux
!! - sfx : salt flux
!! - fr_i : ice fraction
!! - tn_ice : sea-ice surface temperature
!! - alb_ice : sea-ice albedo (recomputed only for coupled mode)
!!
!! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
!! Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
!! These refs are now obsolete since everything has been revised
!! The ref should be Rousset et al., 2015
!!---------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! number of iteration
!
INTEGER :: ji, jj, jl, jk ! dummy loop indices
REAL(wp) :: zqsr ! New solar flux received by the ocean
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d ! 2D workspace
!!---------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('iceupdate')
IF( kt == nit000 .AND. lwp ) THEN
WRITE(numout,*)
WRITE(numout,*)'ice_update_flx: update fluxes (mass, salt and heat) at the ice-ocean interface'
WRITE(numout,*)'~~~~~~~~~~~~~~'
ENDIF
! Net heat flux on top of the ice-ocean (W.m-2)
!----------------------------------------------
IF( ln_cndflx ) THEN ! ice-atm interface = conduction (and melting) fluxes
DO_2D( 0, 0, 0, 0 )
qt_atm_oi(ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * ( qns_oce(ji,jj) + qsr_oce(ji,jj) ) + qemp_oce(ji,jj) &
& + SUM( a_i_b(ji,jj,:) * ( qcn_ice(ji,jj,:) + qml_ice(ji,jj,:) + qtr_ice_top(ji,jj,:) ) ) &
& + qemp_ice(ji,jj)
END_2D
DO_2D( 0, 0, 0, 0 )
qt_atm_oi(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj)
END_2D
ENDIF
! --- case we bypass ice thermodynamics --- !
IF( .NOT. ln_icethd ) THEN ! we suppose ice is impermeable => ocean is isolated from atmosphere
DO_2D( 0, 0, 0, 0 )
qt_atm_oi (ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * ( qns_oce(ji,jj) + qsr_oce(ji,jj) ) + qemp_oce(ji,jj)
qt_oce_ai (ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj)
emp_ice (ji,jj) = 0._wp
qemp_ice (ji,jj) = 0._wp
qevap_ice (ji,jj,:) = 0._wp
END_2D
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
! Solar heat flux reaching the ocean (max) = zqsr (W.m-2)
!---------------------------------------------------
IF( ln_cndflx ) THEN ! ice-atm interface = conduction (and melting) fluxes
zqsr = ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) + SUM( a_i_b (ji,jj,:) * qtr_ice_bot(ji,jj,:) )
ELSE ! ice-atm interface = solar and non-solar fluxes
zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - qtr_ice_bot(ji,jj,:) ) )
ENDIF
! Total heat flux reaching the ocean = qt_oce_ai (W.m-2)
!---------------------------------------------------
IF( ln_icethd ) THEN
qt_oce_ai(ji,jj) = qt_atm_oi(ji,jj) - hfx_sum(ji,jj) - hfx_bom(ji,jj) - hfx_bog(ji,jj) &
& - hfx_dif(ji,jj) - hfx_opw(ji,jj) - hfx_snw(ji,jj) &
& + hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) &
& + hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) + hfx_spr(ji,jj)
ENDIF
! New qsr and qns used to compute the oceanic heat flux at the next time step
!----------------------------------------------------------------------------
! if warming and some ice remains, then we suppose that the whole solar flux has been consumed to melt the ice
! else ( cooling or no ice left ), then we suppose that no solar flux has been consumed
!
IF( fhld(ji,jj) > 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN !-- warming and some ice remains
! solar flux transmitted thru the 1st level of the ocean (i.e. not used by sea-ice)
qsr(ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * ( 1._wp - frq_m(ji,jj) ) &
! + solar flux transmitted thru ice and the 1st ocean level (also not used by sea-ice)
& + SUM( a_i_b(ji,jj,:) * qtr_ice_bot(ji,jj,:) ) * ( 1._wp - frq_m(ji,jj) )
!
ELSE !-- cooling or no ice left
qsr(ji,jj) = zqsr
ENDIF
!
! the non-solar is simply derived from the solar flux
qns(ji,jj) = qt_oce_ai(ji,jj) - qsr(ji,jj)
! Mass flux at the atm. surface
!-----------------------------------
wfx_sub(ji,jj) = wfx_snw_sub(ji,jj) + wfx_ice_sub(ji,jj)
! Mass flux at the ocean surface
!------------------------------------
! ice-ocean mass flux
wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) &
& + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj)
! snw-ocean mass flux
wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj)
! total mass flux at the ocean/ice interface
fmmflx(ji,jj) = - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_pnd(ji,jj) - wfx_err_sub(ji,jj) ! ice-ocean mass flux saved at least for biogeochemical model
emp (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_pnd(ji,jj) - wfx_err_sub(ji,jj) ! atm-ocean + ice-ocean mass flux
! Salt flux at the ocean surface
!------------------------------------------
sfx(ji,jj) = sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj) &
& + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj)
! Mass of snow and ice per unit area
!----------------------------------------
snwice_mass_b(ji,jj) = snwice_mass(ji,jj) ! save mass from the previous ice time step
! ! new mass per unit area
snwice_mass (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) + rhow * (vt_ip(ji,jj) + vt_il(ji,jj)) )
END_2D
CALL lbc_lnk( 'iceupdate', snwice_mass, 'T', 1.0_wp, snwice_mass_b, 'T', 1.0_wp ) ! needed for sshwzv and dynspg_ts (lbc on emp is done in sbcmod)
! time evolution of snow+ice mass
snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_Dt_ice
! Storing the transmitted variables
!----------------------------------
fr_i (:,:) = at_i(:,:) ! Sea-ice fraction
tn_ice(:,:,:) = t_su(A2D(0),:) ! Ice surface temperature
! Snow/ice albedo (only if sent to coupler, useless in forced mode)
!------------------------------------------------------------------
CALL ice_alb( ln_pnd_alb, t_su(A2D(0),:), h_i(A2D(0),:), h_s(A2D(0),:), a_ip_eff(:,:,:), h_ip(A2D(0),:), cloud_fra(:,:), & ! <<== in
& alb_ice(:,:,:) ) ! ==>> out
!
IF( lrst_ice ) THEN !* write snwice_mass fields in the restart file
CALL update_rst( 'WRITE', kt )
ENDIF
!
! output all fluxes
!------------------
!
! --- salt fluxes [kg/m2/s] --- !
! ! sfxice = sfxbog + sfxbom + sfxsum + sfxsni + sfxopw + sfxres + sfxdyn + sfxbri + sfxsub + sfxlam
IF( iom_use('sfxice' ) ) CALL iom_put( 'sfxice', sfx * 1.e-03 ) ! salt flux from total ice growth/melt
IF( iom_use('sfxbog' ) ) CALL iom_put( 'sfxbog', sfx_bog * 1.e-03 ) ! salt flux from bottom growth
IF( iom_use('sfxbom' ) ) CALL iom_put( 'sfxbom', sfx_bom * 1.e-03 ) ! salt flux from bottom melting
IF( iom_use('sfxsum' ) ) CALL iom_put( 'sfxsum', sfx_sum * 1.e-03 ) ! salt flux from surface melting
IF( iom_use('sfxlam' ) ) CALL iom_put( 'sfxlam', sfx_lam * 1.e-03 ) ! salt flux from lateral melting
IF( iom_use('sfxsni' ) ) CALL iom_put( 'sfxsni', sfx_sni * 1.e-03 ) ! salt flux from snow ice formation
IF( iom_use('sfxopw' ) ) CALL iom_put( 'sfxopw', sfx_opw * 1.e-03 ) ! salt flux from open water formation
IF( iom_use('sfxdyn' ) ) CALL iom_put( 'sfxdyn', sfx_dyn * 1.e-03 ) ! salt flux from ridging rafting
IF( iom_use('sfxbri' ) ) CALL iom_put( 'sfxbri', sfx_bri * 1.e-03 ) ! salt flux from brines
IF( iom_use('sfxsub' ) ) CALL iom_put( 'sfxsub', sfx_sub * 1.e-03 ) ! salt flux from sublimation
IF( iom_use('sfxres' ) ) CALL iom_put( 'sfxres', sfx_res(A2D(0)) * 1.e-03 ) ! salt flux from undiagnosed processes
! --- mass fluxes [kg/m2/s] --- !
CALL iom_put( 'emp_oce', emp_oce ) ! emp over ocean (taking into account the snow blown away from the ice)
CALL iom_put( 'emp_ice', emp_ice ) ! emp over ice (taking into account the snow blown away from the ice)
! ! vfxice = vfxbog + vfxbom + vfxsum + vfxsni + vfxopw + vfxdyn + vfxres + vfxlam + vfxpnd
CALL iom_put( 'vfxice' , wfx_ice ) ! mass flux from total ice growth/melt
CALL iom_put( 'vfxbog' , wfx_bog ) ! mass flux from bottom growth
CALL iom_put( 'vfxbom' , wfx_bom ) ! mass flux from bottom melt
CALL iom_put( 'vfxsum' , wfx_sum ) ! mass flux from surface melt
CALL iom_put( 'vfxlam' , wfx_lam ) ! mass flux from lateral melt
CALL iom_put( 'vfxsni' , wfx_sni ) ! mass flux from snow-ice formation
CALL iom_put( 'vfxopw' , wfx_opw ) ! mass flux from growth in open water
CALL iom_put( 'vfxdyn' , wfx_dyn ) ! mass flux from dynamics (ridging)
CALL iom_put( 'vfxpnd' , wfx_pnd ) ! mass flux from melt ponds
CALL iom_put( 'vfxsub' , wfx_ice_sub ) ! mass flux from ice sublimation (ice-atm.)
CALL iom_put( 'vfxsub_err', wfx_err_sub ) ! "excess" of sublimation sent to ocean
CALL iom_put( 'vfxres' , wfx_res(A2D(0)) ) ! mass flux from undiagnosed processes
IF ( iom_use( 'vfxthin' ) ) THEN ! mass flux from ice growth in open water + thin ice (<20cm) => comparable to observations
WHERE( hm_i(:,:) < 0.2 .AND. hm_i(:,:) > 0. ) ; z2d = wfx_bog
ELSEWHERE ; z2d = 0._wp
END WHERE
CALL iom_put( 'vfxthin', wfx_opw + z2d )
DEALLOCATE( z2d )
ENDIF
! ! vfxsnw = vfxsnw_sni + vfxsnw_dyn + vfxsnw_sum
CALL iom_put( 'vfxsnw' , wfx_snw ) ! mass flux from total snow growth/melt
CALL iom_put( 'vfxsnw_sum' , wfx_snw_sum ) ! mass flux from snow melt at the surface
CALL iom_put( 'vfxsnw_sni' , wfx_snw_sni ) ! mass flux from snow melt during snow-ice formation
CALL iom_put( 'vfxsnw_dyn' , wfx_snw_dyn ) ! mass flux from dynamics (ridging)
CALL iom_put( 'vfxsnw_sub' , wfx_snw_sub ) ! mass flux from snow sublimation (ice-atm.)
CALL iom_put( 'vfxsnw_pre' , wfx_spr ) ! snow precip
! --- heat fluxes [W/m2] --- !
! ! qt_atm_oi - qt_oce_ai = hfxdhc - ( dihctrp + dshctrp )
IF( iom_use('qsr_oce' ) ) CALL iom_put( 'qsr_oce' , qsr_oce * ( 1._wp - at_i_b ) ) ! solar flux at ocean surface
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
IF( iom_use('qsr_ice' ) ) CALL iom_put( 'qsr_ice' , SUM( qsr_ice * a_i_b, dim=3 ) ) ! solar flux at ice surface
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
IF( iom_use('qtr_ice_bot') ) CALL iom_put( 'qtr_ice_bot', SUM( qtr_ice_bot * a_i_b, dim=3 ) ) ! solar flux transmitted thru ice
IF( iom_use('qtr_ice_top') ) CALL iom_put( 'qtr_ice_top', SUM( qtr_ice_top * a_i_b, dim=3 ) ) ! solar flux transmitted thru ice surface
IF( iom_use('qt_oce' ) ) CALL iom_put( 'qt_oce' , ( qsr_oce + qns_oce ) * ( 1._wp - at_i_b ) + qemp_oce )
IF( iom_use('qt_ice' ) ) CALL iom_put( 'qt_ice' , SUM( ( qns_ice + qsr_ice ) * a_i_b, dim=3 ) + qemp_ice )
IF( iom_use('qt_oce_ai' ) ) CALL iom_put( 'qt_oce_ai' , qt_oce_ai * smask0 ) ! total heat flux at the ocean surface: interface oce-(ice+atm)
IF( iom_use('qt_atm_oi' ) ) CALL iom_put( 'qt_atm_oi' , qt_atm_oi * smask0 ) ! total heat flux at the oce-ice surface: interface atm-(ice+oce)
IF( iom_use('qemp_oce' ) ) CALL iom_put( 'qemp_oce' , qemp_oce ) ! Downward Heat Flux from E-P over ocean
IF( iom_use('qemp_ice' ) ) CALL iom_put( 'qemp_ice' , qemp_ice ) ! Downward Heat Flux from E-P over ice
! heat fluxes from ice transformations
! ! hfxdhc = hfxbog + hfxbom + hfxsum + hfxopw + hfxdif + hfxsnw - ( hfxthd + hfxdyn + hfxres + hfxsub + hfxspr )
CALL iom_put ('hfxbog' , hfx_bog ) ! heat flux used for ice bottom growth
CALL iom_put ('hfxbom' , hfx_bom ) ! heat flux used for ice bottom melt
CALL iom_put ('hfxsum' , hfx_sum ) ! heat flux used for ice surface melt
CALL iom_put ('hfxopw' , hfx_opw ) ! heat flux used for ice formation in open water
CALL iom_put ('hfxdif' , hfx_dif ) ! heat flux used for ice temperature change
CALL iom_put ('hfxsnw' , hfx_snw ) ! heat flux used for snow melt
CALL iom_put ('hfxerr' , hfx_err_dif ) ! heat flux error after heat diffusion
! heat fluxes associated with mass exchange (freeze/melt/precip...)
CALL iom_put ('hfxthd' , hfx_thd ) !
CALL iom_put ('hfxdyn' , hfx_dyn ) !
CALL iom_put ('hfxsub' , hfx_sub ) !
CALL iom_put ('hfxspr' , hfx_spr ) ! Heat flux from snow precip heat content
CALL iom_put ('hfxres' , hfx_res(A2D(0)) ) !
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
! other heat fluxes
IF( iom_use('hfxsensib' ) ) CALL iom_put( 'hfxsensib' , qsb_ice_bot * at_i_b ) ! Sensible oceanic heat flux
IF( iom_use('hfxcndbot' ) ) CALL iom_put( 'hfxcndbot' , SUM( qcn_ice_bot * a_i_b, dim=3 ) ) ! Bottom conduction flux
IF( iom_use('hfxcndtop' ) ) CALL iom_put( 'hfxcndtop' , SUM( qcn_ice_top * a_i_b, dim=3 ) ) ! Surface conduction flux
IF( iom_use('hfxmelt' ) ) CALL iom_put( 'hfxmelt' , SUM( qml_ice * a_i_b, dim=3 ) ) ! Surface melt flux
IF( iom_use('hfxldmelt' ) ) CALL iom_put( 'hfxldmelt' , fhld * at_i_b ) ! Heat in lead for ice melting
IF( iom_use('hfxldgrow' ) ) CALL iom_put( 'hfxldgrow' , qlead * r1_Dt_ice ) ! Heat in lead for ice growth
! controls
!---------
#if ! defined key_agrif
IF( ln_icediachk ) CALL ice_cons_final('iceupdate') ! conservation
#endif
IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints
IF( sn_cfctl%l_prtctl ) CALL ice_prt3D ('iceupdate') ! prints
IF( ln_timing ) CALL timing_stop ('iceupdate') ! timing
!
END SUBROUTINE ice_update_flx
SUBROUTINE ice_update_tau( kt, pu_oce, pv_oce )
!!-------------------------------------------------------------------
!! *** ROUTINE ice_update_tau ***
!!
!! ** Purpose : Update the ocean surface stresses due to the ice
!!
!! ** Action : * at each ice time step (every nn_fsbc time step):
!! - compute the modulus of ice-ocean relative velocity
!! tmod_io = rhoco * | U_ice-U_oce |
!! - update the modulus of stress at ocean surface
!! taum = (1-a) * taum + a * tmod_io * | U_ice-U_oce |
!! * at each ocean time step (every kt):
!! compute linearized ice-ocean stresses as
!! Utau = tmod_io * | U_ice - pU_oce |
!! using instantaneous current ocean velocity (usually before)
!!
!! NB: - ice-ocean rotation angle no more allowed
!! - here we make an approximation: taum is only computed every ice time step
!! This avoids mutiple average to pass from U,V grids to T grids
!! taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
!! ** Outputs : - utau, vtau : surface ocean i- and j-stress (T-pts) updated with ice-ocean fluxes
!! - taum : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
!!---------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! ocean time-step index
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pu_oce, pv_oce ! surface ocean currents
!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zutau_ice, zu_t, zmodt ! local scalar
REAL(wp) :: zvtau_ice, zv_t, zrhoco ! - -
REAL(wp) :: zflagi ! - -
!!---------------------------------------------------------------------
Clement Rousset
committed
IF( ln_timing ) CALL timing_start('iceupdate')
IF( kt == nit000 .AND. lwp ) THEN
WRITE(numout,*)
WRITE(numout,*)'ice_update_tau: update stress at the ice-ocean interface'
WRITE(numout,*)'~~~~~~~~~~~~~~'
ENDIF
zrhoco = rho0 * rn_cio
!
IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN !== Ice time-step only ==! (i.e. surface module time-step)
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* rhoco * |U_ice-U_oce| at T-point
! ! 2*(U_ice-U_oce) at T-point
zu_t = ( u_ice(ji,jj) + u_ice(ji-1,jj) ) - ( u_oce(ji,jj) + u_oce(ji-1,jj) ) ! u_oce = ssu_m add () for
zv_t = ( v_ice(ji,jj) + v_ice(ji,jj-1) ) - ( v_oce(ji,jj) + v_oce(ji,jj-1) ) ! v_oce = ssv_m NP repro
! ! |U_ice-U_oce|^2
zmodt = 0.25_wp * ( zu_t * zu_t + zv_t * zv_t )
!
tmod_io(ji,jj) = zrhoco * SQRT( zmodt )
END_2D
IF( nn_hls == 1 ) CALL lbc_lnk( 'iceupdate', tmod_io, 'T', 1._wp )
!
DO_2D( 0, 0, 0, 0 ) !* save the air-ocean stresses at ice time-step
zu_t = ( u_ice(ji,jj) + u_ice(ji-1,jj) ) - ( u_oce(ji,jj) + u_oce(ji-1,jj) ) ! u_oce = ssu_m add () for
zv_t = ( v_ice(ji,jj) + v_ice(ji,jj-1) ) - ( v_oce(ji,jj) + v_oce(ji,jj-1) ) ! v_oce = ssv_m NP repro
zmodt = 0.25_wp * ( zu_t * zu_t + zv_t * zv_t )
! ! update the ocean stress modulus
taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * zrhoco * zmodt
!
utau_oce(ji,jj) = utau(ji,jj)
vtau_oce(ji,jj) = vtau(ji,jj)
END_2D
!
ENDIF
!
! !== every ocean time-step ==!
IF ( ln_drgice_imp ) THEN
! Save drag with right sign to update top drag in the ocean implicit friction
DO_2D( 1, 1, 1, 1 )
rCdU_ice(ji,jj) = -r1_rho0 * tmod_io(ji,jj) * at_i(ji,jj) * tmask(ji,jj,1)
END_2D
zflagi = 0._wp
ELSE
zflagi = 1._wp
ENDIF
!
DO_2D( 0, 0, 0, 0 ) !* update the stress WITHOUT an ice-ocean rotation angle
! ! linearized quadratic drag formulation
zutau_ice = 0.5_wp * tmod_io(ji,jj) &
& * ( ( u_ice(ji,jj) + u_ice(ji-1,jj) ) - ( pu_oce(ji,jj) + pu_oce(ji-1,jj) ) ) ! add () for
zvtau_ice = 0.5_wp * tmod_io(ji,jj) &
& * ( ( v_ice(ji,jj) + v_ice(ji,jj-1) ) - ( pv_oce(ji,jj) + pv_oce(ji,jj-1) ) ) ! NP repro
utau(ji,jj) = ( 1._wp - at_i(ji,jj) ) * utau_oce(ji,jj) + at_i(ji,jj) * zutau_ice
vtau(ji,jj) = ( 1._wp - at_i(ji,jj) ) * vtau_oce(ji,jj) + at_i(ji,jj) * zvtau_ice
Clement Rousset
committed
IF( ln_timing ) CALL timing_stop('iceupdate')
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
!
END SUBROUTINE ice_update_tau
SUBROUTINE ice_update_init
!!-------------------------------------------------------------------
!! *** ROUTINE ice_update_init ***
!!
!! ** Purpose : allocate ice-ocean stress fields and read restarts
!! containing the snow & ice mass
!!
!!-------------------------------------------------------------------
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zcoefu, zcoefv, zcoeff ! local scalar
!!-------------------------------------------------------------------
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'ice_update_init: ice-ocean stress init'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
!
! ! allocate ice_update array
IF( ice_update_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'ice_update_init : unable to allocate standard arrays' )
!
CALL update_rst( 'READ' ) !* read or initialize all required files
!
END SUBROUTINE ice_update_init
SUBROUTINE update_rst( cdrw, kt )
!!---------------------------------------------------------------------
!! *** ROUTINE rhg_evp_rst ***
!!
!! ** Purpose : Read or write RHG file in restart file
!!
!! ** Method : use of IOM library
!!----------------------------------------------------------------------
CHARACTER(len=*) , INTENT(in) :: cdrw ! 'READ'/'WRITE' flag
INTEGER, OPTIONAL, INTENT(in) :: kt ! ice time-step
!
INTEGER :: iter ! local integer
INTEGER :: id1 ! local integer
!!----------------------------------------------------------------------
!
IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialize
! ! ---------------
IF( ln_rstart ) THEN !* Read the restart file
!
id1 = iom_varid( numrir, 'snwice_mass' , ldstop = .FALSE. )
!
IF( id1 > 0 ) THEN ! fields exist
CALL iom_get( numrir, jpdom_auto, 'snwice_mass' , snwice_mass )
CALL iom_get( numrir, jpdom_auto, 'snwice_mass_b', snwice_mass_b )
ELSE ! start from rest
IF(lwp) WRITE(numout,*) ' ==>> previous run without snow-ice mass output then set it'
snwice_mass (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) + rhow * (vt_ip(:,:) + vt_il(:,:)) )
snwice_mass_b(:,:) = snwice_mass(:,:)
ENDIF
ELSE !* Start from rest
IF(lwp) WRITE(numout,*) ' ==>> start from rest: set the snow-ice mass'
snwice_mass (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) + rhow * (vt_ip(:,:) + vt_il(:,:)) )
snwice_mass_b(:,:) = snwice_mass(:,:)
ENDIF
!
ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
! ! -------------------
IF(lwp) WRITE(numout,*) '---- update-rst ----'
iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 1
!
CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass' , snwice_mass )
CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b )
!
ENDIF
!
END SUBROUTINE update_rst
#else
!!----------------------------------------------------------------------
!! Default option Dummy module NO SI3 sea-ice model
!!----------------------------------------------------------------------
#endif
!!======================================================================
END MODULE iceupdate