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
58
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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
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
195
196
197
198
199
200
MODULE icedyn_adv_umx
!!==============================================================================
!! *** MODULE icedyn_adv_umx ***
!! sea-ice : advection using the ULTIMATE-MACHO scheme
!!==============================================================================
!! History : 3.6 ! 2014-11 (C. Rousset, G. Madec) Original code
!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
!!----------------------------------------------------------------------
#if defined key_si3
!!----------------------------------------------------------------------
!! 'key_si3' SI3 sea-ice model
!!----------------------------------------------------------------------
!! ice_dyn_adv_umx : update the tracer fields
!! ultimate_x(_y) : compute a tracer value at velocity points using ULTIMATE scheme at various orders
!! macho : compute the fluxes
!! nonosc_ice : limit the fluxes using a non-oscillatory algorithm
!!----------------------------------------------------------------------
USE phycst ! physical constant
USE dom_oce ! ocean domain
USE sbc_oce , ONLY : nn_fsbc ! update frequency of surface boundary condition
USE ice ! sea-ice variables
USE icevar ! sea-ice: operations
!
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)
IMPLICIT NONE
PRIVATE
PUBLIC ice_dyn_adv_umx ! called by icedyn_adv.F90
!
INTEGER, PARAMETER :: np_advS = 1 ! advection for S and T: dVS/dt = -div( uVS ) => np_advS = 1
! or dVS/dt = -div( uA * uHS / u ) => np_advS = 2
! or dVS/dt = -div( uV * uS / u ) => np_advS = 3
INTEGER, PARAMETER :: np_limiter = 1 ! limiter: 1 = nonosc
! 2 = superbee
! 3 = h3
LOGICAL :: ll_upsxy = .TRUE. ! alternate directions for upstream
LOGICAL :: ll_hoxy = .TRUE. ! alternate directions for high order
LOGICAL :: ll_neg = .TRUE. ! if T interpolated at u/v points is negative or v_i < 1.e-6
! then interpolate T at u/v points using the upstream scheme
LOGICAL :: ll_prelim = .FALSE. ! prelimiter from: Zalesak(1979) eq. 14 => not well defined in 2D
!
REAL(wp) :: z1_6 = 1._wp / 6._wp ! =1/6
REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120
!
INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: imsk_small, jmsk_small
!
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: icedyn_adv_umx.F90 15049 2021-06-23 16:17:30Z clem $
!! Software governed by the CeCILL licence (./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, &
& pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
!!----------------------------------------------------------------------
!! *** ROUTINE ice_dyn_adv_umx ***
!!
!! ** Purpose : Compute the now trend due to total advection of
!! tracers and add it to the general trend of tracer equations
!! using an "Ultimate-Macho" scheme
!!
!! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)
INTEGER , INTENT(in ) :: kt ! time step
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity
REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_i ! ice thickness
REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_s ! snw thickness
REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_ip ! ice pond thickness
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond concentration
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_il ! melt pond lid volume
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
!
INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices
INTEGER :: icycle ! number of sub-timestep for the advection
REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers
REAL(wp) :: zdt, z1_dt, zvi_cen
REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication
REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box
REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2
REAL(wp), DIMENSION(jpi,jpj,jpl) :: zu_cat, zv_cat
REAL(wp), DIMENSION(jpi,jpj,jpl) :: zua_ho, zva_ho, zua_ups, zva_ups
REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai , z1_aip, zhvar
REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max, zs_i, zsi_max
REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: ze_i, zei_max
REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: ze_s, zes_max
!
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs
!! diagnostics
REAL(wp), DIMENSION(jpi,jpj) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat
!!----------------------------------------------------------------------
!
IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme'
!
! --- Record max of the surrounding 9-pts (for call Hbig) --- !
! thickness and salinity
WHERE( pv_i(:,:,:) >= epsi10 ) ; zs_i(:,:,:) = psv_i(:,:,:) / pv_i(:,:,:)
ELSEWHERE ; zs_i(:,:,:) = 0._wp
END WHERE
CALL icemax3D( ph_i , zhi_max )
CALL icemax3D( ph_s , zhs_max )
CALL icemax3D( ph_ip, zhip_max)
CALL icemax3D( zs_i , zsi_max )
CALL lbc_lnk( 'icedyn_adv_umx', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp )
!
! enthalpies
DO jk = 1, nlay_i
WHERE( pv_i(:,:,:) >= epsi10 ) ; ze_i(:,:,jk,:) = pe_i(:,:,jk,:) / pv_i(:,:,:)
ELSEWHERE ; ze_i(:,:,jk,:) = 0._wp
END WHERE
END DO
DO jk = 1, nlay_s
WHERE( pv_s(:,:,:) >= epsi10 ) ; ze_s(:,:,jk,:) = pe_s(:,:,jk,:) / pv_s(:,:,:)
ELSEWHERE ; ze_s(:,:,jk,:) = 0._wp
END WHERE
END DO
CALL icemax4D( ze_i , zei_max )
CALL icemax4D( ze_s , zes_max )
CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_umx', zes_max, 'T', 1._wp )
!
!
! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- !
! Note: the advection split is applied at the next time-step in order to avoid blocking global comm.
! this should not affect too much the stability
zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rDt_ice * r1_e1u(:,:) )
zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rDt_ice * r1_e2v(:,:) ) )
! non-blocking global communication send zcflnow and receive zcflprv
CALL mpp_delay_max( 'icedyn_adv_umx', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
IF( zcflprv(1) > .5 ) THEN ; icycle = 2
ELSE ; icycle = 1
ENDIF
zdt = rDt_ice / REAL(icycle)
z1_dt = 1._wp / zdt
! --- transport --- !
zudy(:,:) = pu_ice(:,:) * e2u(:,:)
zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
!
! setup transport for each ice cat
DO jl = 1, jpl
zu_cat(:,:,jl) = zudy(:,:)
zv_cat(:,:,jl) = zvdx(:,:)
END DO
!
! --- define velocity for advection: u*grad(H) --- !
DO_2D( nn_hls-1, nn_hls, nn_hls, nn_hls )
IF ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN ; zcu_box(ji,jj) = 0._wp
ELSEIF( pu_ice(ji,jj) > 0._wp ) THEN ; zcu_box(ji,jj) = pu_ice(ji-1,jj)
ELSE ; zcu_box(ji,jj) = pu_ice(ji ,jj)
ENDIF
END_2D
DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls )
IF ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN ; zcv_box(ji,jj) = 0._wp
ELSEIF( pv_ice(ji,jj) > 0._wp ) THEN ; zcv_box(ji,jj) = pv_ice(ji,jj-1)
ELSE ; zcv_box(ji,jj) = pv_ice(ji,jj )
ENDIF
END_2D
!---------------!
!== advection ==!
!---------------!
DO jt = 1, icycle
! diagnostics
zdiag_adv_mass(:,:) = SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos &
& + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow
zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi
zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 )
! record at_i before advection (for open water)
zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
! inverse of A and Ap
WHERE( pa_i(:,:,:) >= epsi20 ) ; z1_ai(:,:,:) = 1._wp / pa_i(:,:,:)
ELSEWHERE ; z1_ai(:,:,:) = 0.
END WHERE
WHERE( pa_ip(:,:,:) >= epsi20 ) ; z1_aip(:,:,:) = 1._wp / pa_ip(:,:,:)
ELSEWHERE ; z1_aip(:,:,:) = 0.
END WHERE
Loading
Loading full blame...