Changeset 5965 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
- Timestamp:
- 2015-12-01T16:35:30+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r4333 r5965 22 22 USE thd_ice ! LIM thermodynamics 23 23 USE dom_ice ! LIM domain 24 USE par_ice ! LIM parameters25 24 USE ice ! LIM variables 26 25 USE limtab ! LIM 2D <==> 1D … … 29 28 USE lib_mpp ! MPP library 30 29 USE wrk_nemo ! work arrays 30 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 31 31 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 32 USE limthd_ent 33 USE limvar 32 34 33 35 IMPLICIT NONE … … 35 37 36 38 PUBLIC lim_thd_lac ! called by lim_thd 37 38 REAL(wp) :: epsi10 = 1.e-10_wp !39 REAL(wp) :: zzero = 0._wp !40 REAL(wp) :: zone = 1._wp !41 39 42 40 !!---------------------------------------------------------------------- … … 71 69 !! - Computation of variation of ice volume and mass 72 70 !! - Computation of frldb after lateral accretion and 73 !! update ht_s_ b, ht_i_band tbif_1d(:,:)71 !! update ht_s_1d, ht_i_1d and tbif_1d(:,:) 74 72 !!------------------------------------------------------------------------ 75 INTEGER :: ji,jj,jk,jl ,jm! dummy loop indices76 INTEGER :: layer, nbpac! local integers77 INTEGER :: ii, ij, iter ! - -78 REAL(wp) :: ztmelts, zdv, z qold, zfrazb, zweight, zalphai, zindb, zinda, zde! local scalars73 INTEGER :: ji,jj,jk,jl ! dummy loop indices 74 INTEGER :: nbpac ! local integers 75 INTEGER :: ii, ij, iter ! - - 76 REAL(wp) :: ztmelts, zdv, zfrazb, zweight, zde ! local scalars 79 77 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - - 80 78 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - - 81 79 LOGICAL :: iterate_frazil ! iterate frazil ice collection thickness 82 80 CHARACTER (len = 15) :: fieldid 83 ! 84 INTEGER , POINTER, DIMENSION(:) :: zcatac ! indexes of categories where new ice grows 81 82 REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean) 83 REAL(wp) :: zEi ! sea ice specific enthalpy (J/kg) 84 REAL(wp) :: zEw ! seawater specific enthalpy (J/kg) 85 REAL(wp) :: zfmdt ! mass flux x time step (kg/m2, >0 towards ocean) 86 87 REAL(wp) :: zv_newfra 88 89 INTEGER , POINTER, DIMENSION(:) :: jcat ! indexes of categories where new ice grows 85 90 REAL(wp), POINTER, DIMENSION(:) :: zswinew ! switch for new ice or not 86 91 … … 93 98 REAL(wp), POINTER, DIMENSION(:) :: zdv_res ! residual volume in case of excessive heat budget 94 99 REAL(wp), POINTER, DIMENSION(:) :: zda_res ! residual area in case of excessive heat budget 95 REAL(wp), POINTER, DIMENSION(:) :: zat_i_ac ! total ice fraction 96 REAL(wp), POINTER, DIMENSION(:) :: zat_i_lev ! total ice fraction for level ice only (type 1) 97 REAL(wp), POINTER, DIMENSION(:) :: zdh_frazb ! accretion of frazil ice at the ice bottom 98 REAL(wp), POINTER, DIMENSION(:) :: zvrel_ac ! relative ice / frazil velocity (1D vector) 99 100 REAL(wp), POINTER, DIMENSION(:,:) :: zhice_old ! previous ice thickness 101 REAL(wp), POINTER, DIMENSION(:,:) :: zdummy ! dummy thickness of new ice 102 REAL(wp), POINTER, DIMENSION(:,:) :: zdhicbot ! thickness of new ice which is accreted vertically 103 REAL(wp), POINTER, DIMENSION(:,:) :: zv_old ! old volume of ice in category jl 104 REAL(wp), POINTER, DIMENSION(:,:) :: za_old ! old area of ice in category jl 105 REAL(wp), POINTER, DIMENSION(:,:) :: za_i_ac ! 1-D version of a_i 106 REAL(wp), POINTER, DIMENSION(:,:) :: zv_i_ac ! 1-D version of v_i 107 REAL(wp), POINTER, DIMENSION(:,:) :: zoa_i_ac ! 1-D version of oa_i 108 REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_ac ! 1-D version of smv_i 109 110 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_ac !: 1-D version of e_i 111 112 REAL(wp), POINTER, DIMENSION(:) :: zqbgow ! heat budget of the open water (negative) 113 REAL(wp), POINTER, DIMENSION(:) :: zdhex ! excessively thick accreted sea ice (hlead-hice) 114 115 REAL(wp), POINTER, DIMENSION(:,:,:) :: zqm0 ! old layer-system heat content 116 REAL(wp), POINTER, DIMENSION(:,:,:) :: zthick0 ! old ice thickness 117 118 REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories 119 REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories 120 REAL(wp), POINTER, DIMENSION(:,:) :: et_i_init, et_i_final ! ice energy summed over categories 121 REAL(wp), POINTER, DIMENSION(:,:) :: et_s_init ! snow energy summed over categories 100 REAL(wp), POINTER, DIMENSION(:) :: zat_i_1d ! total ice fraction 101 REAL(wp), POINTER, DIMENSION(:) :: zv_frazb ! accretion of frazil ice at the ice bottom 102 REAL(wp), POINTER, DIMENSION(:) :: zvrel_1d ! relative ice / frazil velocity (1D vector) 103 104 REAL(wp), POINTER, DIMENSION(:,:) :: zv_b ! old volume of ice in category jl 105 REAL(wp), POINTER, DIMENSION(:,:) :: za_b ! old area of ice in category jl 106 REAL(wp), POINTER, DIMENSION(:,:) :: za_i_1d ! 1-D version of a_i 107 REAL(wp), POINTER, DIMENSION(:,:) :: zv_i_1d ! 1-D version of v_i 108 REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_1d ! 1-D version of smv_i 109 110 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_1d !: 1-D version of e_i 111 122 112 REAL(wp), POINTER, DIMENSION(:,:) :: zvrel ! relative ice / frazil velocity 113 114 REAL(wp) :: zcai = 1.4e-3_wp 123 115 !!-----------------------------------------------------------------------! 124 116 125 CALL wrk_alloc( jpij, zcatac) ! integer117 CALL wrk_alloc( jpij, jcat ) ! integer 126 118 CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 127 CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex ) 128 CALL wrk_alloc( jpij,jpl, zhice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac ) 129 CALL wrk_alloc( jpij,jkmax,jpl, ze_i_ac ) 130 CALL wrk_alloc( jpij,jkmax+1,jpl, zqm0, zthick0 ) 131 CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel ) 132 133 et_i_init(:,:) = 0._wp 134 et_s_init(:,:) = 0._wp 135 vt_i_init(:,:) = 0._wp 136 vt_s_init(:,:) = 0._wp 137 138 !------------------------------------------------------------------------------! 139 ! 1) Conservation check and changes in each ice category 140 !------------------------------------------------------------------------------! 141 IF( con_i ) THEN 142 CALL lim_column_sum ( jpl, v_i , vt_i_init) 143 CALL lim_column_sum ( jpl, v_s , vt_s_init) 144 CALL lim_column_sum_energy ( jpl, nlay_i , e_i , et_i_init) 145 CALL lim_column_sum ( jpl, e_s(:,:,1,:) , et_s_init) 146 ENDIF 147 119 CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 120 CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d ) 121 CALL wrk_alloc( jpij,nlay_i,jpl, ze_i_1d ) 122 CALL wrk_alloc( jpi,jpj, zvrel ) 123 124 CALL lim_var_agg(1) 125 CALL lim_var_glo2eqv 148 126 !------------------------------------------------------------------------------| 149 127 ! 2) Convert units for ice internal energy … … 154 132 DO ji = 1, jpi 155 133 !Energy of melting q(S,T) [J.m-3] 156 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * REAL( nlay_i ) 157 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 158 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 134 rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 ) ) !0 if no ice 135 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl), epsi20 ) * REAL( nlay_i, wp ) 159 136 END DO 160 137 END DO … … 179 156 180 157 ! Default new ice thickness 181 hicol(:,:) = hiccrit(1)182 183 IF( fraz_swi == 1._wp) THEN158 hicol(:,:) = rn_hnewice 159 160 IF( ln_frazil ) THEN 184 161 185 162 !-------------------- … … 190 167 zhicrit = 0.04 ! frazil ice thickness 191 168 ztwogp = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav 192 zsqcd = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag)169 zsqcd = 1.0 / SQRT( 1.3 * zcai ) ! 1/SQRT(airdensity*drag) 193 170 zgamafr = 0.03 194 171 195 DO jj = 1, jpj 196 DO ji = 1, jpi 197 198 IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN 172 DO jj = 2, jpj 173 DO ji = 2, jpi 174 IF ( qlead(ji,jj) < 0._wp ) THEN 199 175 !------------- 200 176 ! Wind stress 201 177 !------------- 202 178 ! C-grid wind stress components 203 ztaux = ( utau_ice(ji-1,jj ) * tmu(ji-1,jj) &204 & + utau_ice(ji ,jj ) * tmu(ji ,jj) ) * 0.5_wp205 ztauy = ( vtau_ice(ji ,jj-1) * tmv(ji ,jj-1) &206 & + vtau_ice(ji ,jj ) * tmv(ji ,jj) ) * 0.5_wp179 ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) & 180 & + utau_ice(ji ,jj ) * umask(ji ,jj ,1) ) * 0.5_wp 181 ztauy = ( vtau_ice(ji ,jj-1) * vmask(ji ,jj-1,1) & 182 & + vtau_ice(ji ,jj ) * vmask(ji ,jj ,1) ) * 0.5_wp 207 183 ! Square root of wind stress 208 ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy) )184 ztenagm = SQRT( SQRT( ztaux**2 + ztauy**2 ) ) 209 185 210 186 !--------------------- 211 187 ! Frazil ice velocity 212 188 !--------------------- 213 zvfrx = zgamafr * zsqcd * ztaux / MAX(ztenagm,epsi10) 214 zvfry = zgamafr * zsqcd * ztauy / MAX(ztenagm,epsi10) 189 rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 190 zvfrx = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 191 zvfry = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 215 192 216 193 !------------------- … … 218 195 !------------------- 219 196 ! C-grid ice velocity 220 zindb = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) ) ) 221 zvgx = zindb * ( u_ice(ji-1,jj ) * tmu(ji-1,jj ) & 222 & + u_ice(ji,jj ) * tmu(ji ,jj ) ) * 0.5_wp 223 zvgy = zindb * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1) & 224 & + v_ice(ji,jj ) * tmv(ji ,jj ) ) * 0.5_wp 197 rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) ) ) 198 zvgx = rswitch * ( u_ice(ji-1,jj ) * umask(ji-1,jj ,1) + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 199 zvgy = rswitch * ( v_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 225 200 226 201 !----------------------------------- … … 248 223 iterate_frazil = .true. 249 224 250 DO WHILE ( iter .LT.100 .AND. iterate_frazil )225 DO WHILE ( iter < 100 .AND. iterate_frazil ) 251 226 zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) & 252 227 - hicol(ji,jj) * zhicrit * ztwogp * zvrel2 … … 264 239 END DO ! loop on ji ends 265 240 END DO ! loop on jj ends 241 ! 242 CALL lbc_lnk( zvrel(:,:), 'T', 1. ) 243 CALL lbc_lnk( hicol(:,:), 'T', 1. ) 266 244 267 245 ENDIF ! End of computation of frazil ice collection thickness … … 276 254 ! This occurs if open water energy budget is negative 277 255 nbpac = 0 256 npac(:) = 0 257 ! 278 258 DO jj = 1, jpj 279 259 DO ji = 1, jpi 280 IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) >0._wp ) THEN260 IF ( qlead(ji,jj) < 0._wp ) THEN 281 261 nbpac = nbpac + 1 282 262 npac( nbpac ) = (jj - 1) * jpi + ji … … 287 267 ! debug point to follow 288 268 jiindex_1d = 0 289 IF( ln_ nicep) THEN290 DO ji = mi0( jiindx), mi1(jiindx)291 DO jj = mj0(j jindx), mj1(jjindx)292 IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) >0._wp ) THEN269 IF( ln_icectl ) THEN 270 DO ji = mi0(iiceprt), mi1(iiceprt) 271 DO jj = mj0(jiceprt), mj1(jiceprt) 272 IF ( qlead(ji,jj) < 0._wp ) THEN 293 273 jiindex_1d = (jj - 1) * jpi + ji 294 274 ENDIF … … 297 277 ENDIF 298 278 299 IF( ln_ nicep) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac279 IF( ln_icectl ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 300 280 301 281 !------------------------------ … … 307 287 IF ( nbpac > 0 ) THEN 308 288 309 CALL tab_2d_1d( nbpac, zat_i_ ac(1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) )289 CALL tab_2d_1d( nbpac, zat_i_1d (1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) ) 310 290 DO jl = 1, jpl 311 CALL tab_2d_1d( nbpac, za_i_ac (1:nbpac,jl), a_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 312 CALL tab_2d_1d( nbpac, zv_i_ac (1:nbpac,jl), v_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 313 CALL tab_2d_1d( nbpac, zoa_i_ac (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 314 CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 291 CALL tab_2d_1d( nbpac, za_i_1d (1:nbpac,jl), a_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 292 CALL tab_2d_1d( nbpac, zv_i_1d (1:nbpac,jl), v_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 293 CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 315 294 DO jk = 1, nlay_i 316 CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 317 END DO ! jk 318 END DO ! jl 319 320 CALL tab_2d_1d( nbpac, qldif_1d (1:nbpac) , qldif , jpi, jpj, npac(1:nbpac) ) 321 CALL tab_2d_1d( nbpac, qcmif_1d (1:nbpac) , qcmif , jpi, jpj, npac(1:nbpac) ) 322 CALL tab_2d_1d( nbpac, t_bo_b (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) ) 323 CALL tab_2d_1d( nbpac, sfx_thd_1d(1:nbpac) , sfx_thd, jpi, jpj, npac(1:nbpac) ) 324 CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac) , rdm_ice, jpi, jpj, npac(1:nbpac) ) 325 CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 326 CALL tab_2d_1d( nbpac, zvrel_ac (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 295 CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 296 END DO 297 END DO 298 299 CALL tab_2d_1d( nbpac, qlead_1d (1:nbpac) , qlead , jpi, jpj, npac(1:nbpac) ) 300 CALL tab_2d_1d( nbpac, t_bo_1d (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) ) 301 CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac) , sfx_opw, jpi, jpj, npac(1:nbpac) ) 302 CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw, jpi, jpj, npac(1:nbpac) ) 303 CALL tab_2d_1d( nbpac, hicol_1d (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 304 CALL tab_2d_1d( nbpac, zvrel_1d (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 305 306 CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac) , hfx_thd, jpi, jpj, npac(1:nbpac) ) 307 CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac) , hfx_opw, jpi, jpj, npac(1:nbpac) ) 327 308 328 309 !------------------------------------------------------------------------------! … … 330 311 !------------------------------------------------------------------------------! 331 312 313 !----------------------------------------- 314 ! Keep old ice areas and volume in memory 315 !----------------------------------------- 316 zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:) 317 za_b(1:nbpac,:) = za_i_1d(1:nbpac,:) 332 318 !---------------------- 333 319 ! Thickness of new ice 334 320 !---------------------- 335 321 DO ji = 1, nbpac 336 zh_newice(ji) = hiccrit(1)337 END DO 338 IF( fraz_swi == 1.0 ) zh_newice(:) = hicol_b(:)322 zh_newice(ji) = rn_hnewice 323 END DO 324 IF( ln_frazil ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 339 325 340 326 !---------------------- 341 327 ! Salinity of new ice 342 328 !---------------------- 343 344 SELECT CASE ( num_sal ) 329 SELECT CASE ( nn_icesal ) 345 330 CASE ( 1 ) ! Sice = constant 346 zs_newice( :) = bulk_sal331 zs_newice(1:nbpac) = rn_icesal 347 332 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 348 333 DO ji = 1, nbpac 349 334 ii = MOD( npac(ji) - 1 , jpi ) + 1 350 335 ij = ( npac(ji) - 1 ) / jpi + 1 351 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij) )336 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_m(ii,ij) ) 352 337 END DO 353 338 CASE ( 3 ) ! Sice = F(z) [multiyear ice] 354 zs_newice( :) = 2.3339 zs_newice(1:nbpac) = 2.3 355 340 END SELECT 356 357 341 358 342 !------------------------- … … 361 345 ! We assume that new ice is formed at the seawater freezing point 362 346 DO ji = 1, nbpac 363 ztmelts = - tmut * zs_newice(ji) + rtt ! Melting point (K) 364 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) & 365 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) ) & 366 & - rcp * ( ztmelts - rtt ) ) 367 ze_newice(ji) = MAX( ze_newice(ji) , 0._wp ) & 368 & + MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) ) * rhoic * lfus 369 END DO ! ji 347 ztmelts = - tmut * zs_newice(ji) + rt0 ! Melting point (K) 348 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_1d(ji) ) & 349 & + lfus * ( 1.0 - ( ztmelts - rt0 ) / MIN( t_bo_1d(ji) - rt0, -epsi10 ) ) & 350 & - rcp * ( ztmelts - rt0 ) ) 351 END DO 352 370 353 !---------------- 371 354 ! Age of new ice … … 373 356 DO ji = 1, nbpac 374 357 zo_newice(ji) = 0._wp 375 END DO ! ji 376 377 !-------------------------- 378 ! Open water energy budget 379 !-------------------------- 380 DO ji = 1, nbpac 381 zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji) !<0 382 END DO ! ji 358 END DO 383 359 384 360 !------------------- … … 386 362 !------------------- 387 363 DO ji = 1, nbpac 388 zv_newice(ji) = - zqbgow(ji) / ze_newice(ji) 364 365 zEi = - ze_newice(ji) * r1_rhoic ! specific enthalpy of forming ice [J/kg] 366 367 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_1d [J/kg] 368 ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 369 370 zdE = zEi - zEw ! specific enthalpy difference [J/kg] 371 372 zfmdt = - qlead_1d(ji) / zdE ! Fm.dt [kg/m2] (<0) 373 ! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point 374 zv_newice(ji) = - zfmdt * r1_rhoic 375 376 zQm = zfmdt * zEw ! heat to the ocean >0 associated with mass flux 377 378 ! Contribution to heat flux to the ocean [W.m-2], >0 379 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice 380 ! Total heat flux used in this process [W.m-2] 381 hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice 382 ! mass flux 383 wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_rdtice 384 ! salt flux 385 sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice 389 386 390 387 ! A fraction zfrazb of frazil ice is accreted at the ice bottom 391 zfrazb = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 392 zdh_frazb(ji) = zfrazb * zv_newice(ji) 388 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 389 zfrazb = rswitch * ( TANH ( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb 390 zv_frazb(ji) = zfrazb * zv_newice(ji) 393 391 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 394 392 END DO 395 396 !------------------------------------397 ! Diags for energy conservation test398 !------------------------------------399 DO ji = 1, nbpac400 ii = MOD( npac(ji) - 1 , jpi ) + 1401 ij = ( npac(ji) - 1 ) / jpi + 1402 !403 zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji)404 !405 vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji) ! volume406 et_i_init(ii,ij) = et_i_init(ii,ij) + zde ! Energy407 408 END DO409 410 ! keep new ice volume in memory411 CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , jpi, jpj )412 393 413 394 !----------------- … … 415 396 !----------------- 416 397 DO ji = 1, nbpac 417 ii = MOD( npac(ji) - 1 , jpi ) + 1418 ij = ( npac(ji) - 1 ) / jpi + 1419 398 za_newice(ji) = zv_newice(ji) / zh_newice(ji) 420 diag_lat_gr(ii,ij) = diag_lat_gr(ii,ij) + zv_newice(ji) * r1_rdtice ! clem 421 END DO !ji 399 END DO 422 400 423 401 !------------------------------------------------------------------------------! … … 425 403 !------------------------------------------------------------------------------! 426 404 427 !----------------------------------------- 428 ! Keep old ice areas and volume in memory 429 !----------------------------------------- 430 zv_old(:,:) = zv_i_ac(:,:) 431 za_old(:,:) = za_i_ac(:,:) 432 433 !------------------------------------------- 434 ! Compute excessive new ice area and volume 435 !------------------------------------------- 405 !------------------------ 406 ! 6.1) lateral ice growth 407 !------------------------ 436 408 ! If lateral ice growth gives an ice concentration gt 1, then 437 409 ! we keep the excessive volume in memory and attribute it later to bottom accretion 438 410 DO ji = 1, nbpac 439 IF ( za_newice(ji) > ( amax - zat_i_ac(ji) ) ) THEN440 zda_res(ji) = za_newice(ji) - ( amax - zat_i_ac(ji) )411 IF ( za_newice(ji) > ( rn_amax - zat_i_1d(ji) ) ) THEN 412 zda_res(ji) = za_newice(ji) - ( rn_amax - zat_i_1d(ji) ) 441 413 zdv_res(ji) = zda_res (ji) * zh_newice(ji) 442 414 za_newice(ji) = za_newice(ji) - zda_res (ji) … … 446 418 zdv_res(ji) = 0._wp 447 419 ENDIF 448 END DO ! ji 449 450 !------------------------------------------------ 451 ! Laterally redistribute new ice volume and area 452 !------------------------------------------------ 453 zat_i_ac(:) = 0._wp 420 END DO 421 422 ! find which category to fill 423 zat_i_1d(:) = 0._wp 454 424 DO jl = 1, jpl 455 425 DO ji = 1, nbpac 456 IF( hi_max (jl-1) < zh_newice(ji) .AND. & 457 & zh_newice(ji) <= hi_max (jl) ) THEN 458 za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji) 459 zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji) 460 zat_i_ac(ji) = zat_i_ac(ji) + za_i_ac (ji,jl) 461 zcatac (ji) = jl 426 IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN 427 za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji) 428 zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji) 429 jcat (ji) = jl 462 430 ENDIF 463 END DO 464 END DO 465 466 !---------------------------------- 467 ! Heat content - lateral accretion 468 !---------------------------------- 469 DO ji = 1, nbpac 470 jl = zcatac(ji) ! categroy in which new ice is put 471 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) ) ! zindb=1 if ice =0 otherwise 472 zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb ! old ice thickness 473 zdhex (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) ) ! difference in thickness 474 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) ) ! ice totally new in jl category 431 zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d (ji,jl) 432 END DO 433 END DO 434 435 ! Heat content 436 DO ji = 1, nbpac 437 jl = jcat(ji) ! categroy in which new ice is put 438 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) ) ! 0 if old ice 475 439 END DO 476 440 477 441 DO jk = 1, nlay_i 478 442 DO ji = 1, nbpac 479 jl = zcatac(ji) 480 zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ] 481 zalphai = MIN( zhice_old(ji,jl) * REAL( jk ) / REAL( nlay_i ), zh_newice(ji) ) & 482 & - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) ) 483 ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji) & 484 + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / REAL( nlay_i ) & 485 + za_newice(ji) * ze_newice(ji) * zalphai & 486 + za_newice(ji) * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) ) 487 END DO 488 END DO 489 490 !----------------------------------------------- 491 ! Add excessive volume of new ice at the bottom 492 !----------------------------------------------- 493 ! If the ice concentration exceeds 1, the remaining volume of new ice 494 ! is equally redistributed among all ice categories in which there is 495 ! ice 496 497 ! Fraction of level ice 498 jm = 1 499 zat_i_lev(:) = 0._wp 500 501 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 502 DO ji = 1, nbpac 503 zat_i_lev(ji) = zat_i_lev(ji) + za_i_ac(ji,jl) 504 END DO 505 END DO 506 507 IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl) 508 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 509 DO ji = 1, nbpac 510 zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) ) 511 zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi10 ) ) ! clem 512 zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi10 ) 513 END DO 514 END DO 515 IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl) 516 517 !--------------------------------- 518 ! Heat content - bottom accretion 519 !--------------------------------- 520 jm = 1 521 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 522 DO ji = 1, nbpac 523 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) ) ! zindb=1 if ice =0 otherwise 524 zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb 525 zdhicbot (ji,jl) = zdv_res(ji) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb & 526 & + zindb * zdh_frazb(ji) ! frazil ice may coalesce 527 zdummy(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb ! thickness of residual ice 528 END DO 529 END DO 530 531 ! old layers thicknesses and enthalpies 532 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 443 jl = jcat(ji) 444 rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 445 ze_i_1d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + & 446 & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) ) & 447 & * rswitch / MAX( zv_i_1d(ji,jl), epsi20 ) 448 END DO 449 END DO 450 451 !------------------------------------------------ 452 ! 6.2) bottom ice growth + ice enthalpy remapping 453 !------------------------------------------------ 454 DO jl = 1, jpl 455 456 ! for remapping 457 h_i_old (1:nbpac,0:nlay_i+1) = 0._wp 458 qh_i_old(1:nbpac,0:nlay_i+1) = 0._wp 533 459 DO jk = 1, nlay_i 534 460 DO ji = 1, nbpac 535 zthick0(ji,jk,jl) = zhice_old(ji,jl) / REAL( nlay_i )536 zqm0 (ji,jk,jl) = ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl)461 h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i 462 qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk) 537 463 END DO 538 464 END DO 539 END DO 540 !!gm ??? why the previous do loop if ocerwriten by the following one ? 541 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 465 466 ! new volumes including lateral/bottom accretion + residual 542 467 DO ji = 1, nbpac 543 zthick0(ji,nlay_i+1,jl) = zdhicbot(ji,jl) 544 zqm0 (ji,nlay_i+1,jl) = ze_newice(ji) * zdhicbot(ji,jl) 545 END DO ! ji 546 END DO ! jl 547 548 ! Redistributing energy on the new grid 549 ze_i_ac(:,:,:) = 0._wp 550 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 551 DO jk = 1, nlay_i 552 DO layer = 1, nlay_i + 1 553 DO ji = 1, nbpac 554 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) 555 ! Redistributing energy on the new grid 556 zweight = MAX ( MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) ) & 557 & - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp ) & 558 & /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb 559 ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl) 560 END DO ! ji 561 END DO ! layer 562 END DO ! jk 563 END DO ! jl 564 565 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 566 DO jk = 1, nlay_i 567 DO ji = 1, nbpac 568 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) 569 ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) & 570 & / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb 571 END DO 572 END DO 573 END DO 574 575 !------------ 576 ! Update age 577 !------------ 578 DO jl = 1, jpl 579 DO ji = 1, nbpac 580 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes 581 zoa_i_ac(ji,jl) = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb 582 END DO 583 END DO 468 rswitch = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) ) 469 zv_newfra = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 ) 470 za_i_1d(ji,jl) = rswitch * za_i_1d(ji,jl) 471 zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra 472 ! for remapping 473 h_i_old (ji,nlay_i+1) = zv_newfra 474 qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra 475 ENDDO 476 ! --- Ice enthalpy remapping --- ! 477 CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 478 ENDDO 584 479 585 480 !----------------- 586 481 ! Update salinity 587 482 !----------------- 588 !clem IF( num_sal == 2 ) THEN589 DO jl = 1, jpl590 DO ji = 1, nbpac591 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes592 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl)593 zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif594 END DO595 END DO596 !clem ENDIF597 598 !--------------------------------599 ! Update mass/salt fluxes (clem)600 !--------------------------------601 483 DO jl = 1, jpl 602 484 DO ji = 1, nbpac 603 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes 604 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl) 605 rdm_ice_1d(ji) = rdm_ice_1d(ji) + zdv * rhoic * zindb 606 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zdv * rhoic * zs_newice(ji) * r1_rdtice * zindb 607 END DO 485 zdv = zv_i_1d(ji,jl) - zv_b(ji,jl) 486 zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 487 END DO 608 488 END DO 609 489 610 490 !------------------------------------------------------------------------------! 611 ! 8) Change 2D vectors to 1D vectors491 ! 7) Change 2D vectors to 1D vectors 612 492 !------------------------------------------------------------------------------! 613 493 DO jl = 1, jpl 614 CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_ac (1:nbpac,jl), jpi, jpj ) 615 CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj ) 616 CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj ) 617 !clem IF ( num_sal == 2 ) & 618 CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj ) 494 CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj ) 495 CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj ) 496 CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj ) 619 497 DO jk = 1, nlay_i 620 CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_ac(1:nbpac,jk,jl), jpi, jpj ) 621 END DO 622 END DO 623 CALL tab_1d_2d( nbpac, sfx_thd, npac(1:nbpac), sfx_thd_1d(1:nbpac), jpi, jpj ) 624 CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj ) 498 CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_1d(1:nbpac,jk,jl), jpi, jpj ) 499 END DO 500 END DO 501 CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj ) 502 CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj ) 503 504 CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj ) 505 CALL tab_1d_2d( nbpac, hfx_opw, npac(1:nbpac), hfx_opw_1d(1:nbpac), jpi, jpj ) 625 506 ! 626 507 ENDIF ! nbpac > 0 627 508 628 509 !------------------------------------------------------------------------------! 629 ! 9) Change units for e_i510 ! 8) Change units for e_i 630 511 !------------------------------------------------------------------------------! 631 512 DO jl = 1, jpl 632 DO jk = 1, nlay_i ! heat content in 10^9 Joules 633 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i ) / unit_fac 513 DO jk = 1, nlay_i 514 DO jj = 1, jpj 515 DO ji = 1, jpi 516 ! heat content in J/m2 517 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 518 END DO 519 END DO 634 520 END DO 635 521 END DO 636 522 637 !------------------------------------------------------------------------------|638 ! 10) Conservation check and changes in each ice category639 !------------------------------------------------------------------------------|640 IF( con_i ) THEN641 CALL lim_column_sum (jpl, v_i, vt_i_final)642 fieldid = 'v_i, limthd_lac'643 CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)644 !645 CALL lim_column_sum_energy(jpl, nlay_i, e_i, et_i_final)646 fieldid = 'e_i, limthd_lac'647 CALL lim_cons_check (et_i_final, et_i_final, 1.0e-3, fieldid)648 !649 CALL lim_column_sum (jpl, v_s, vt_s_final)650 fieldid = 'v_s, limthd_lac'651 CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)652 !653 ! CALL lim_column_sum (jpl, e_s(:,:,1,:) , et_s_init)654 ! fieldid = 'e_s, limthd_lac'655 ! CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)656 IF( ln_nicep ) THEN657 DO ji = mi0(jiindx), mi1(jiindx)658 DO jj = mj0(jjindx), mj1(jjindx)659 WRITE(numout,*) ' vt_i_init : ', vt_i_init (ji,jj)660 WRITE(numout,*) ' vt_i_final: ', vt_i_final(ji,jj)661 WRITE(numout,*) ' et_i_init : ', et_i_init (ji,jj)662 WRITE(numout,*) ' et_i_final: ', et_i_final(ji,jj)663 END DO664 END DO665 ENDIF666 !667 ENDIF668 523 ! 669 CALL wrk_dealloc( jpij, zcatac) ! integer524 CALL wrk_dealloc( jpij, jcat ) ! integer 670 525 CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 671 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex ) 672 CALL wrk_dealloc( jpij,jpl, zhice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac ) 673 CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_ac ) 674 CALL wrk_dealloc( jpij,jkmax+1,jpl, zqm0, zthick0 ) 675 CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel ) 526 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 527 CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d ) 528 CALL wrk_dealloc( jpij,nlay_i,jpl, ze_i_1d ) 529 CALL wrk_dealloc( jpi,jpj, zvrel ) 676 530 ! 677 531 END SUBROUTINE lim_thd_lac
Note: See TracChangeset
for help on using the changeset viewer.