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

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_da.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

update trunk with OpenMP parallelization

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