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.
limthd_da.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_da.F90 @ 8369

Last change on this file since 8369 was 8369, checked in by clem, 7 years ago

STEP4 (4): put all thermodynamics in 1D (limitd_th OK)
File size: 12.5 KB
Line 
1MODULE limthd_da
2   !!======================================================================
3   !!                       ***  MODULE limthd_da ***
4   !! LIM-3 sea-ice :  computation of lateral melting in the ice
5   !!======================================================================
6   !! History :   4.0   ! 2016-03 (C. Rousset) original code
7   !!---------------------------------------------------------------------
8#if defined key_lim3
9   !!----------------------------------------------------------------------
10   !!   'key_lim3'                                      LIM-3 sea-ice model
11   !!----------------------------------------------------------------------
12   !!   lim_thd_da   : sea ice lateral melting
13   !!----------------------------------------------------------------------
14   USE par_oce        ! ocean parameters
15   USE phycst         ! physical constants (ocean directory)
16   USE sbc_oce , ONLY : sst_m
17   USE ice            ! LIM variables
18   USE thd_ice        ! thermodynamic sea-ice variables
19   USE limtab         ! 1D <==> 2D transformation
20   !
21   USE lib_mpp        ! MPP library
22   USE wrk_nemo       ! work arrays
23   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   lim_thd_da        ! called by limthd module
29
30   !!----------------------------------------------------------------------
31   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
32   !! $Id: limthd_da.F90 5123 2015-03-04 16:06:03Z clem $
33   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
34   !!----------------------------------------------------------------------
35CONTAINS
36
37   SUBROUTINE lim_thd_da
38      !!-------------------------------------------------------------------
39      !!                ***  ROUTINE lim_thd_da  ***   
40      !!   
41      !! ** Purpose :   computes sea ice lateral melting
42      !!
43      !! ** Method  :   dA/dt = - P * W   [s-1]
44      !!                   W = melting velocity [m.s-1]
45      !!                   P = perimeter of ice-ocean lateral interface normalized by grid cell area [m.m-2]
46      !!
47      !!                   W = m1 * (Tw -Tf)**m2                    --- originally from Josberger 1979 ---
48      !!                      (Tw - Tf) = elevation of water temp above freezing
49      !!                      m1 and m2 = (1.6e-6 , 1.36) best fit from field experiment near the coast of Prince Patrick Island (Perovich 1983) => static ice
50      !!                      m1 and m2 = (3.0e-6 , 1.36) best fit from MIZEX 84 experiment (Maykut and Perovich 1987) => moving ice
51      !!
52      !!                   P = N * pi * D                           --- from Rothrock and Thorndike 1984 ---
53      !!                      D = mean floe caliper diameter
54      !!                      N = number of floes = ice area / floe area(average) = A / (Cs * D**2)
55      !!                         A = ice concentration
56      !!                         Cs = deviation from a square (square:Cs=1 ; circle:Cs=pi/4 ; floe:Cs=0.66)
57      !!
58      !!                   D = Dmin * ( Astar / (Astar-A) )**beta   --- from Lupkes et al., 2012 (eq. 26-27) ---
59      !!                                                             
60      !!                      Astar = 1 / ( 1 - (Dmin/Dmax)**(1/beta) )
61      !!                      Dmin = minimum floe diameter (recommended to be 8m +- 20%)
62      !!                      Dmax = maximum floe diameter (recommended to be 300m, but it does not impact melting much except for Dmax<100m)
63      !!                      beta = 1.0 +-20% (recommended value)
64      !!                           = 0.3 best fit for western Fram Strait and Antarctica
65      !!                           = 1.4 best fit for eastern Fram Strait
66      !!
67      !! ** Tunable parameters  :   We propose to tune the lateral melting via 2 parameters
68      !!                               Dmin [6-10m]   => 6  vs 8m = +40% melting at the peak (A~0.5)
69      !!                                                 10 vs 8m = -20% melting
70      !!                               beta [0.8-1.2] => decrease = more melt and melt peaks toward higher concentration
71      !!                                                                  (A~0.5 for beta=1 ; A~0.8 for beta=0.2)
72      !!                                                 0.3 = best fit for western Fram Strait and Antarctica
73      !!                                                 1.4 = best fit for eastern Fram Strait
74      !!
75      !! ** Note   :   Former and more simple formulations for floe diameters can be found in Mai (1995),
76      !!               Birnbaum and Lupkes (2002), Lupkes and Birnbaum (2005). They are reviewed in Lupkes et al 2012
77      !!               A simpler implementation for CICE can be found in Bitz et al (2001) and Tsamados et al (2015)
78      !!
79      !! ** References
80      !!    Bitz, C. M., Holland, M. M., Weaver, A. J., & Eby, M. (2001).
81      !!              Simulating the ice‐thickness distribution in a coupled climate model.
82      !!              Journal of Geophysical Research: Oceans, 106(C2), 2441-2463.
83      !!    Josberger, E. G. (1979).
84      !!              Laminar and turbulent boundary layers adjacent to melting vertical ice walls in salt water
85      !!              (No. SCIENTIFIC-16). WASHINGTON UNIV SEATTLE DEPT OF ATMOSPHERIC SCIENCES.
86      !!    Lüpkes, C., Gryanik, V. M., Hartmann, J., & Andreas, E. L. (2012).
87      !!              A parametrization, based on sea ice morphology, of the neutral atmospheric drag coefficients
88      !!              for weather prediction and climate models.
89      !!              Journal of Geophysical Research: Atmospheres, 117(D13).
90      !!    Maykut, G. A., & Perovich, D. K. (1987).
91      !!              The role of shortwave radiation in the summer decay of a sea ice cover.
92      !!              Journal of Geophysical Research: Oceans, 92(C7), 7032-7044.
93      !!    Perovich, D. K. (1983).
94      !!              On the summer decay of a sea ice cover. (Doctoral dissertation, University of Washington).
95      !!    Rothrock, D. A., & Thorndike, A. S. (1984).
96      !!              Measuring the sea ice floe size distribution.
97      !!              Journal of Geophysical Research: Oceans, 89(C4), 6477-6486.
98      !!    Tsamados, M., Feltham, D., Petty, A., Schroeder, D., & Flocco, D. (2015).
99      !!              Processes controlling surface, bottom and lateral melt of Arctic sea ice in a state of the art sea ice model.
100      !!              Phil. Trans. R. Soc. A, 373(2052), 20140167.
101      !!---------------------------------------------------------------------
102      INTEGER             ::   ji, jj, jk, jl     ! dummy loop indices
103      REAL(wp) ::   ztmelts             ! local scalar
104      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg)
105      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg)
106      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
107      REAL(wp)            ::   zastar, zdfloe, zperi, zwlat, zda
108      REAL(wp), PARAMETER ::   zdmax = 300._wp
109      REAL(wp), PARAMETER ::   zcs   = 0.66_wp
110      REAL(wp), PARAMETER ::   zm1   = 3.e-6_wp
111      REAL(wp), PARAMETER ::   zm2   = 1.36_wp
112      !
113      REAL(wp), DIMENSION(jpij) ::   zda_tot
114      !!---------------------------------------------------------------------
115
116      ! select ice covered grid points
117      nidx = 0 ; idxice(:) = 0
118      DO jj = 1, jpj
119         DO ji = 1, jpi
120            IF ( at_i(ji,jj) > epsi10 ) THEN
121               nidx         = nidx  + 1
122               idxice(nidx) = (jj - 1) * jpi + ji
123            ENDIF
124         END DO
125      END DO
126
127      IF( nidx > 0 ) THEN
128         !------------------------------------------------------------!
129         ! --- Calculate reduction of total sea ice concentration --- !
130         !------------------------------------------------------------!
131         CALL tab_2d_1d( nidx, idxice(1:nidx), at_i_1d(1:nidx), at_i  )
132         CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d(1:nidx), t_bo  )
133         CALL tab_2d_1d( nidx, idxice(1:nidx), sst_1d (1:nidx), sst_m )
134         
135         zastar = 1._wp / ( 1._wp - (rn_dmin / zdmax)**(1._wp/rn_beta) )
136         DO ji = 1, nidx   
137            zdfloe = rn_dmin * ( zastar / ( zastar - at_i_1d(ji) ) )**rn_beta         ! Mean floe caliper diameter [m]
138            zperi  = at_i_1d(ji) * rpi / ( zcs * zdfloe )                             ! Mean perimeter of the floe = N*pi*D = (A/cs*D^2)*pi*D [m.m-2]
139            zwlat  = zm1 * ( MAX( 0._wp, sst_1d(ji) - ( t_bo_1d(ji) - rt0 ) ) )**zm2  ! Melt speed rate [m/s]
140           
141            zda_tot(ji) = - MIN( zwlat * zperi * rdt_ice, at_i_1d(ji) )               ! sea ice concentration decrease
142         END DO
143         
144         !---------------------------------------------------------------------------------------------!
145         ! --- Distribute reduction among ice categories and calculate associated ice-ocean fluxes --- !
146         !---------------------------------------------------------------------------------------------!
147         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_lam_1d(1:nidx), sfx_lam )
148         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd )
149         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_lam_1d(1:nidx), wfx_lam )
150         
151         DO jl = 1, jpl
152           
153            CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl)  )
154            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d(1:nidx), ht_i(:,:,jl) )
155            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_s_1d(1:nidx), ht_s(:,:,jl) )
156            CALL tab_2d_1d( nidx, idxice(1:nidx), sm_i_1d(1:nidx), sm_i(:,:,jl) )
157            DO jk = 1, nlay_i
158               CALL tab_2d_1d( nidx, idxice(1:nidx), e_i_1d(1:nidx,jk), e_i(:,:,jk,jl) )
159            END DO
160            DO jk = 1, nlay_s
161               CALL tab_2d_1d( nidx, idxice(1:nidx), e_s_1d(1:nidx,jk), e_s(:,:,jk,jl) )
162            END DO
163           
164            DO ji = 1, nidx
165               IF( a_i_1d(ji) > 0._wp ) THEN
166                  ! decrease of concentration for the category jl
167                  !    each category contributes to melting in proportion to its concentration
168                  zda = zda_tot(ji) * a_i_1d(ji) / at_i_1d(ji)
169                 
170                  ! Contribution to salt flux
171                  sfx_lam_1d(ji) = sfx_lam_1d(ji) - rhoic *  ht_i_1d(ji) * zda * sm_i_1d(ji) * r1_rdtice
172                 
173                  ! Contribution to heat flux into the ocean [W.m-2], (<0) 
174                  hfx_thd_1d(ji) = hfx_thd_1d(ji) + zda_tot(ji) / at_i_1d(ji) * SUM( e_i_1d(ji,1:nlay_i) + e_s_1d(ji,1) ) * r1_rdtice 
175                 
176                  ! Contribution to mass flux
177                  wfx_lam_1d(ji) =  wfx_lam_1d(ji) - zda * r1_rdtice * ( rhoic * ht_i_1d(ji) + rhosn * ht_s_1d(ji) )
178                 
179                  !! adjust e_i ???
180                  e_i_1d(ji,1:nlay_i) = e_i_1d(ji,1:nlay_i) * ( 1._wp + zda_tot(ji) / at_i_1d(ji) )
181                  e_s_1d(ji,1)        = e_s_1d(ji,1)        * ( 1._wp + zda_tot(ji) / at_i_1d(ji) )
182                 
183                  ! new concentration
184                  a_i_1d(ji) = a_i_1d(ji) + zda
185                 
186                  ! ensure that ht_i = 0 where a_i = 0
187                  IF( a_i_1d(ji) == 0._wp ) THEN
188                     ht_i_1d(ji) = 0._wp
189                     ht_s_1d(ji) = 0._wp
190                  ENDIF
191               ENDIF
192            END DO
193           
194            CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d    (1:nidx), a_i (:,:,jl) )
195            CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d   (1:nidx), ht_i(:,:,jl) )
196            CALL tab_1d_2d( nidx, idxice(1:nidx), ht_s_1d   (1:nidx), ht_s(:,:,jl) )
197            DO jk = 1, nlay_s
198               CALL tab_1d_2d( nidx, idxice(1:nidx), e_s_1d(1:nidx,jk), e_s(:,:,jk,jl) )
199            END DO
200            DO jk = 1, nlay_i
201               CALL tab_1d_2d( nidx, idxice(1:nidx), e_i_1d(1:nidx,jk), e_i(:,:,jk,jl) )
202            END DO
203           
204         END DO
205         
206         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_lam_1d(1:nidx), sfx_lam )
207         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd )
208         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_lam_1d(1:nidx), wfx_lam )
209         
210      ENDIF
211      !
212   END SUBROUTINE lim_thd_da
213   
214#else
215   !!----------------------------------------------------------------------
216   !!   Default option         Dummy Module          No LIM-3 sea-ice model
217   !!----------------------------------------------------------------------
218CONTAINS
219   SUBROUTINE lim_thd_da          ! Empty routine
220   END SUBROUTINE lim_thd_da
221#endif
222   !!======================================================================
223END MODULE limthd_da
Note: See TracBrowser for help on using the repository browser.