Changeset 6225 for branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r4688 r6225 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 32 USE limthd_ent 33 USE limvar 33 34 34 35 IMPLICIT NONE … … 36 37 37 38 PUBLIC lim_thd_lac ! called by lim_thd 38 39 REAL(wp) :: epsi10 = 1.e-10_wp !40 REAL(wp) :: epsi20 = 1.e-20_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, zfrazb, zweight, z indb, 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 ! - - … … 89 87 REAL(wp) :: zv_newfra 90 88 91 INTEGER , POINTER, DIMENSION(:) :: jcat ! indexes of categories where new ice grows89 INTEGER , POINTER, DIMENSION(:) :: jcat ! indexes of categories where new ice grows 92 90 REAL(wp), POINTER, DIMENSION(:) :: zswinew ! switch for new ice or not 93 91 … … 101 99 REAL(wp), POINTER, DIMENSION(:) :: zda_res ! residual area in case of excessive heat budget 102 100 REAL(wp), POINTER, DIMENSION(:) :: zat_i_1d ! total ice fraction 103 REAL(wp), POINTER, DIMENSION(:) :: zat_i_lev ! total ice fraction for level ice only (type 1) 104 REAL(wp), POINTER, DIMENSION(:) :: zv_frazb ! accretion of frazil ice at the ice bottom 101 REAL(wp), POINTER, DIMENSION(:) :: zv_frazb ! accretion of frazil ice at the ice bottom 105 102 REAL(wp), POINTER, DIMENSION(:) :: zvrel_1d ! relative ice / frazil velocity (1D vector) 106 103 107 REAL(wp), POINTER, DIMENSION(:,:) :: zv_old ! old volume of ice in category jl 108 REAL(wp), POINTER, DIMENSION(:,:) :: za_old ! old area of ice in category jl 109 REAL(wp), POINTER, DIMENSION(:,:) :: za_i_1d ! 1-D version of a_i 110 REAL(wp), POINTER, DIMENSION(:,:) :: zv_i_1d ! 1-D version of v_i 111 REAL(wp), POINTER, DIMENSION(:,:) :: zoa_i_1d ! 1-D version of oa_i 112 REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_1d ! 1-D version of smv_i 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 113 109 114 110 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_1d !: 1-D version of e_i 115 111 116 112 REAL(wp), POINTER, DIMENSION(:,:) :: zvrel ! relative ice / frazil velocity 113 114 REAL(wp) :: zcai = 1.4e-3_wp 117 115 !!-----------------------------------------------------------------------! 118 116 119 117 CALL wrk_alloc( jpij, jcat ) ! integer 120 118 CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 121 CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, z at_i_lev, zv_frazb, zvrel_1d )122 CALL wrk_alloc( jpij,jpl, zv_ old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )123 CALL wrk_alloc( jpij, jkmax,jpl, ze_i_1d )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 ) 124 122 CALL wrk_alloc( jpi,jpj, zvrel ) 125 123 124 CALL lim_var_agg(1) 125 CALL lim_var_glo2eqv 126 126 !------------------------------------------------------------------------------| 127 127 ! 2) Convert units for ice internal energy … … 132 132 DO ji = 1, jpi 133 133 !Energy of melting q(S,T) [J.m-3] 134 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 135 e_i(ji,jj,jk,jl) = zindb * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i ) 136 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 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 ) 137 136 END DO 138 137 END DO … … 157 156 158 157 ! Default new ice thickness 159 hicol(:,:) = hiccrit160 161 IF( fraz_swi == 1) THEN158 hicol(:,:) = rn_hnewice 159 160 IF( ln_frazil ) THEN 162 161 163 162 !-------------------- … … 168 167 zhicrit = 0.04 ! frazil ice thickness 169 168 ztwogp = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav 170 zsqcd = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag)169 zsqcd = 1.0 / SQRT( 1.3 * zcai ) ! 1/SQRT(airdensity*drag) 171 170 zgamafr = 0.03 172 171 173 DO jj = 1, jpj 174 DO ji = 1, jpi 175 172 DO jj = 2, jpj 173 DO ji = 2, jpi 176 174 IF ( qlead(ji,jj) < 0._wp ) THEN 177 175 !------------- … … 179 177 !------------- 180 178 ! C-grid wind stress components 181 ztaux = ( utau_ice(ji-1,jj ) * tmu(ji-1,jj) &182 & + utau_ice(ji ,jj ) * tmu(ji ,jj) ) * 0.5_wp183 ztauy = ( vtau_ice(ji ,jj-1) * tmv(ji ,jj-1) &184 & + 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 185 183 ! Square root of wind stress 186 184 ztenagm = SQRT( SQRT( ztaux**2 + ztauy**2 ) ) … … 189 187 ! Frazil ice velocity 190 188 !--------------------- 191 zindb= MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )192 zvfrx = zindb* zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )193 zvfry = zindb* 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 ) 194 192 195 193 !------------------- … … 197 195 !------------------- 198 196 ! C-grid ice velocity 199 zindb = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) ) ) 200 zvgx = zindb * ( u_ice(ji-1,jj ) * tmu(ji-1,jj ) & 201 & + u_ice(ji,jj ) * tmu(ji ,jj ) ) * 0.5_wp 202 zvgy = zindb * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1) & 203 & + 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 204 200 205 201 !----------------------------------- … … 227 223 iterate_frazil = .true. 228 224 229 DO WHILE ( iter .LT.100 .AND. iterate_frazil )225 DO WHILE ( iter < 100 .AND. iterate_frazil ) 230 226 zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) & 231 227 - hicol(ji,jj) * zhicrit * ztwogp * zvrel2 … … 243 239 END DO ! loop on ji ends 244 240 END DO ! loop on jj ends 241 ! 242 CALL lbc_lnk( zvrel(:,:), 'T', 1. ) 243 CALL lbc_lnk( hicol(:,:), 'T', 1. ) 245 244 246 245 ENDIF ! End of computation of frazil ice collection thickness … … 255 254 ! This occurs if open water energy budget is negative 256 255 nbpac = 0 256 npac(:) = 0 257 ! 257 258 DO jj = 1, jpj 258 259 DO ji = 1, jpi … … 266 267 ! debug point to follow 267 268 jiindex_1d = 0 268 IF( ln_ nicep) THEN269 DO ji = mi0( jiindx), mi1(jiindx)270 DO jj = mj0(j jindx), mj1(jjindx)269 IF( ln_icectl ) THEN 270 DO ji = mi0(iiceprt), mi1(iiceprt) 271 DO jj = mj0(jiceprt), mj1(jiceprt) 271 272 IF ( qlead(ji,jj) < 0._wp ) THEN 272 273 jiindex_1d = (jj - 1) * jpi + ji … … 276 277 ENDIF 277 278 278 IF( ln_ nicep) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac279 IF( ln_icectl ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 279 280 280 281 !------------------------------ … … 290 291 CALL tab_2d_1d( nbpac, za_i_1d (1:nbpac,jl), a_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 291 292 CALL tab_2d_1d( nbpac, zv_i_1d (1:nbpac,jl), v_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 292 CALL tab_2d_1d( nbpac, zoa_i_1d (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) )293 293 CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 294 294 DO jk = 1, nlay_i 295 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 ! jk297 END DO ! jl296 END DO 297 END DO 298 298 299 299 CALL tab_2d_1d( nbpac, qlead_1d (1:nbpac) , qlead , jpi, jpj, npac(1:nbpac) ) 300 CALL tab_2d_1d( nbpac, t_bo_ b(1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) )300 CALL tab_2d_1d( nbpac, t_bo_1d (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) ) 301 301 CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac) , sfx_opw, jpi, jpj, npac(1:nbpac) ) 302 302 CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw, jpi, jpj, npac(1:nbpac) ) 303 CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw, jpi, jpj, npac(1:nbpac) ) 304 CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 303 CALL tab_2d_1d( nbpac, hicol_1d (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 305 304 CALL tab_2d_1d( nbpac, zvrel_1d (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 306 305 … … 315 314 ! Keep old ice areas and volume in memory 316 315 !----------------------------------------- 317 zv_old(:,:) = zv_i_1d(:,:) 318 za_old(:,:) = za_i_1d(:,:) 319 316 zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:) 317 za_b(1:nbpac,:) = za_i_1d(1:nbpac,:) 320 318 !---------------------- 321 319 ! Thickness of new ice 322 320 !---------------------- 323 321 DO ji = 1, nbpac 324 zh_newice(ji) = hiccrit325 END DO 326 IF( fraz_swi == 1 ) 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) 327 325 328 326 !---------------------- 329 327 ! Salinity of new ice 330 328 !---------------------- 331 SELECT CASE ( n um_sal )329 SELECT CASE ( nn_icesal ) 332 330 CASE ( 1 ) ! Sice = constant 333 zs_newice( :) = bulk_sal331 zs_newice(1:nbpac) = rn_icesal 334 332 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 335 333 DO ji = 1, nbpac 336 334 ii = MOD( npac(ji) - 1 , jpi ) + 1 337 335 ij = ( npac(ji) - 1 ) / jpi + 1 338 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) ) 339 337 END DO 340 338 CASE ( 3 ) ! Sice = F(z) [multiyear ice] 341 zs_newice( :) = 2.3339 zs_newice(1:nbpac) = 2.3 342 340 END SELECT 343 341 … … 347 345 ! We assume that new ice is formed at the seawater freezing point 348 346 DO ji = 1, nbpac 349 ztmelts = - tmut * zs_newice(ji) + rt t! Melting point (K)350 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_ b(ji) ) &351 & + lfus * ( 1.0 - ( ztmelts - rt t ) / MIN( t_bo_b(ji) - rtt, -epsi10 ) ) &352 & - rcp * ( ztmelts - rt t) )353 END DO ! ji347 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 354 352 355 353 !---------------- … … 358 356 DO ji = 1, nbpac 359 357 zo_newice(ji) = 0._wp 360 END DO ! ji358 END DO 361 359 362 360 !------------------- … … 365 363 DO ji = 1, nbpac 366 364 367 zEi = - ze_newice(ji) / rhoic! specific enthalpy of forming ice [J/kg]368 369 zEw = rcp * ( t_bo_ b(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_b[J/kg]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] 370 368 ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 371 369 … … 374 372 zfmdt = - qlead_1d(ji) / zdE ! Fm.dt [kg/m2] (<0) 375 373 ! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point 376 zv_newice(ji) = - zfmdt /rhoic374 zv_newice(ji) = - zfmdt * r1_rhoic 377 375 378 376 zQm = zfmdt * zEw ! heat to the ocean >0 associated with mass flux … … 388 386 389 387 ! A fraction zfrazb of frazil ice is accreted at the ice bottom 390 zinda= 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )391 zfrazb = zinda * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 *maxfrazb388 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 392 390 zv_frazb(ji) = zfrazb * zv_newice(ji) 393 391 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) … … 411 409 ! we keep the excessive volume in memory and attribute it later to bottom accretion 412 410 DO ji = 1, nbpac 413 IF ( za_newice(ji) > ( amax - zat_i_1d(ji) ) ) THEN414 zda_res(ji) = za_newice(ji) - ( amax - zat_i_1d(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) ) 415 413 zdv_res(ji) = zda_res (ji) * zh_newice(ji) 416 414 za_newice(ji) = za_newice(ji) - zda_res (ji) … … 438 436 DO ji = 1, nbpac 439 437 jl = jcat(ji) ! categroy in which new ice is put 440 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_ old(ji,jl) ) ) ! 0 if old ice438 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) ) ! 0 if old ice 441 439 END DO 442 440 … … 444 442 DO ji = 1, nbpac 445 443 jl = jcat(ji) 446 zinda= MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )444 rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 447 445 ze_i_1d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + & 448 & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_ old(ji,jl) ) &449 & * zinda/ MAX( zv_i_1d(ji,jl), epsi20 )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 ) 450 448 END DO 451 449 END DO … … 461 459 DO jk = 1, nlay_i 462 460 DO ji = 1, nbpac 463 h_i_old (ji,jk) = zv_i_1d(ji,jl) / REAL( nlay_i )461 h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i 464 462 qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk) 465 463 END DO … … 468 466 ! new volumes including lateral/bottom accretion + residual 469 467 DO ji = 1, nbpac 470 zinda= MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) )471 zv_newfra = zinda* ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 )472 za_i_1d(ji,jl) = zinda* za_i_1d(ji,jl)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) 473 471 zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra 474 475 472 ! for remapping 476 473 h_i_old (ji,nlay_i+1) = zv_newfra 477 474 qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra 478 475 ENDDO 479 480 476 ! --- Ice enthalpy remapping --- ! 481 IF( zv_newfra > 0._wp ) THEN 482 CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 483 ENDIF 484 477 CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 485 478 ENDDO 486 487 !------------488 ! Update age489 !------------490 DO jl = 1, jpl491 DO ji = 1, nbpac492 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) ) ! 0 if no ice and 1 if yes493 zoa_i_1d(ji,jl) = za_old(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb494 END DO495 END DO496 479 497 480 !----------------- … … 500 483 DO jl = 1, jpl 501 484 DO ji = 1, nbpac 502 zdv = zv_i_1d(ji,jl) - zv_ old(ji,jl)485 zdv = zv_i_1d(ji,jl) - zv_b(ji,jl) 503 486 zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 504 487 END DO … … 511 494 CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj ) 512 495 CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj ) 513 CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_1d(1:nbpac,jl), jpi, jpj )514 496 CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj ) 515 497 DO jk = 1, nlay_i … … 518 500 END DO 519 501 CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj ) 520 CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj )521 502 CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj ) 522 503 … … 533 514 DO jj = 1, jpj 534 515 DO ji = 1, jpi 535 ! heat content in J oules536 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ) * unit_fac )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 537 518 END DO 538 519 END DO … … 543 524 CALL wrk_dealloc( jpij, jcat ) ! integer 544 525 CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 545 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, z at_i_lev, zv_frazb, zvrel_1d )546 CALL wrk_dealloc( jpij,jpl, zv_ old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )547 CALL wrk_dealloc( jpij, jkmax,jpl, ze_i_1d )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 ) 548 529 CALL wrk_dealloc( jpi,jpj, zvrel ) 549 530 !
Note: See TracChangeset
for help on using the changeset viewer.
