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 @ 8342

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

simplify the code
File size: 12.3 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      !------------------------------------------------------------!
128      ! --- Calculate reduction of total sea ice concentration --- !
129      !------------------------------------------------------------!
130      zastar = 1._wp / ( 1._wp - (rn_dmin / zdmax)**(1._wp/rn_beta) )
131
132      CALL tab_2d_1d( nidx, idxice(1:nidx), at_i_1d(1:nidx), at_i  )
133      CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d(1:nidx), t_bo  )
134      CALL tab_2d_1d( nidx, idxice(1:nidx), sst_1d (1:nidx), sst_m )
135
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      DO jl = 1, jpl
148
149         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d    (1:nidx), a_i(:,:,jl)  )
150         CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d   (1:nidx), ht_i(:,:,jl) )
151         CALL tab_2d_1d( nidx, idxice(1:nidx), ht_s_1d   (1:nidx), ht_s(:,:,jl) )
152         CALL tab_2d_1d( nidx, idxice(1:nidx), sm_i_1d   (1:nidx), sm_i(:,:,jl) )
153         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_lam_1d(1:nidx), sfx_lam      )
154         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd      )
155         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_lam_1d(1:nidx), wfx_lam      )
156         DO jk = 1, nlay_i
157            CALL tab_2d_1d( nidx, idxice(1:nidx), e_i_1d(1:nidx,jk), e_i(:,:,jk,jl) )
158         END DO
159         DO jk = 1, nlay_s
160            CALL tab_2d_1d( nidx, idxice(1:nidx), e_s_1d(1:nidx,jk), e_s(:,:,jk,jl) )
161         END DO
162
163         DO ji = 1, nidx
164            IF( a_i_1d(ji) > epsi10 ) THEN
165               ! decrease of concentration for the category jl
166               !    each category contributes to melting in proportion to its concentration
167               zda = zda_tot(ji) * a_i_1d(ji) / at_i_1d(ji)
168               
169               ! Contribution to salt flux
170               sfx_lam_1d(ji) = sfx_lam_1d(ji) - rhoic *  ht_i_1d(ji) * zda * sm_i_1d(ji) * r1_rdtice
171               
172               ! Contribution to heat flux into the ocean [W.m-2], (<0) 
173               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zda / a_i_1d(ji) * SUM( e_i_1d(ji,:) + e_s_1d(ji,1) ) * r1_rdtice 
174             
175               ! Contribution to mass flux
176               wfx_lam_1d(ji) =  wfx_lam_1d(ji) - zda * r1_rdtice * ( rhoic * ht_i_1d(ji) + rhosn * ht_s_1d(ji) )
177           
178               !! adjust e_i ???
179               e_i_1d(ji,:) = e_i_1d(ji,:) * ( 1._wp + zda / a_i_1d(ji) )
180               e_s_1d(ji,1) = e_s_1d(ji,1) * ( 1._wp + zda / a_i_1d(ji) )
181 
182               ! new concentration
183               a_i_1d(ji) = a_i_1d(ji) + zda
184
185               ! ensure that ht_i = 0 where a_i = 0
186               IF( a_i_1d(ji) == 0._wp ) THEN
187                  ht_i_1d(ji) = 0._wp
188                  ht_s_1d(ji) = 0._wp
189               ENDIF
190
191            ENDIF     
192           
193         END DO
194
195!! je pense qu'il faut ajuster e_i mais je ne sais pas comment
196         DO jk = 1, nlay_s
197            CALL tab_1d_2d( nidx, idxice(1:nidx), e_s_1d(1:nidx,jk), e_s(:,:,jk,jl) )
198         END DO
199         DO jk = 1, nlay_i
200            CALL tab_1d_2d( nidx, idxice(1:nidx), e_i_1d(1:nidx,jk), e_i(:,:,jk,jl) )
201         END DO
202         
203         CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d    (1:nidx), a_i (:,:,jl) )
204         CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d   (1:nidx), ht_i(:,:,jl) )
205         CALL tab_1d_2d( nidx, idxice(1:nidx), ht_s_1d   (1:nidx), ht_s(:,:,jl) )
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      END DO
211           
212      !
213   END SUBROUTINE lim_thd_da
214   
215#else
216   !!----------------------------------------------------------------------
217   !!   Default option         Dummy Module          No LIM-3 sea-ice model
218   !!----------------------------------------------------------------------
219CONTAINS
220   SUBROUTINE lim_thd_da          ! Empty routine
221   END SUBROUTINE lim_thd_da
222#endif
223   !!======================================================================
224END MODULE limthd_da
Note: See TracBrowser for help on using the repository browser.