Changeset 4902 for branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3
- Timestamp:
- 2014-11-27T17:13:38+01:00 (10 years ago)
- Location:
- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 22 edited
-
ice.F90 (modified) (12 diffs)
-
iceini.F90 (modified) (4 diffs)
-
limcons.F90 (modified) (2 diffs)
-
limdiahsb.F90 (modified) (1 diff)
-
limdyn.F90 (modified) (4 diffs)
-
limistate.F90 (modified) (1 diff)
-
limitd_me.F90 (modified) (6 diffs)
-
limitd_th.F90 (modified) (19 diffs)
-
limsbc.F90 (modified) (5 diffs)
-
limthd.F90 (modified) (7 diffs)
-
limthd_dh.F90 (modified) (34 diffs)
-
limthd_dif.F90 (modified) (44 diffs)
-
limthd_ent.F90 (modified) (1 diff)
-
limthd_lac.F90 (modified) (15 diffs)
-
limthd_sal.F90 (modified) (5 diffs)
-
limtrp.F90 (modified) (6 diffs)
-
limupdate1.F90 (modified) (6 diffs)
-
limupdate2.F90 (modified) (7 diffs)
-
limvar.F90 (modified) (5 diffs)
-
limwri.F90 (modified) (1 diff)
-
par_ice.F90 (modified) (2 diffs)
-
thd_ice.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r4900 r4902 105 105 !! ** Global variables | 106 106 !!-------------|-------------|---------------------------------|-------| 107 !! a_i | a_i_ b| Ice concentration | |107 !! a_i | a_i_1d | Ice concentration | | 108 108 !! v_i | - | Ice volume per unit area | m | 109 109 !! v_s | - | Snow volume per unit area | m | … … 111 111 !! oa_i ! - ! Sea ice areal age content | day | 112 112 !! e_i ! - ! Ice enthalpy | 10^9 J| 113 !! - ! q_i_ b! Ice enthalpy per unit vol. | J/m3 |113 !! - ! q_i_1d ! Ice enthalpy per unit vol. | J/m3 | 114 114 !! e_s ! - ! Snow enthalpy | 10^9 J| 115 !! - ! q_s_ b! Snow enthalpy per unit vol. | J/m3 |115 !! - ! q_s_1d ! Snow enthalpy per unit vol. | J/m3 | 116 116 !! | 117 117 !!-------------|-------------|---------------------------------|-------| … … 120 120 !!-------------|-------------|---------------------------------|-------| 121 121 !! | 122 !! ht_i | ht_i_ b| Ice thickness | m |123 !! ht_s ! ht_s_ b| Snow depth | m |124 !! sm_i ! sm_i_ b| Sea ice bulk salinity ! ppt |125 !! s_i ! s_i_ b| Sea ice salinity profile ! ppt |122 !! ht_i | ht_i_1d | Ice thickness | m | 123 !! ht_s ! ht_s_1d | Snow depth | m | 124 !! sm_i ! sm_i_1d | Sea ice bulk salinity ! ppt | 125 !! s_i ! s_i_1d | Sea ice salinity profile ! ppt | 126 126 !! o_i ! - | Sea ice Age ! days | 127 !! t_i ! t_i_ b| Sea ice temperature ! K |128 !! t_s ! t_s_ b| Snow temperature ! K |129 !! t_su ! t_su_ b| Sea ice surface temperature ! K |127 !! t_i ! t_i_1d | Sea ice temperature ! K | 128 !! t_s ! t_s_1d | Snow temperature ! K | 129 !! t_su ! t_su_1d | Sea ice surface temperature ! K | 130 130 !! | 131 131 !! notes: the ice model only sees a bulk (i.e., vertically averaged) | … … 142 142 !! *** Category-summed state variables (diagnostic) *** | 143 143 !! ******************************************************************* | 144 !! at_i | at_i_ b| Total ice concentration | |144 !! at_i | at_i_1d | Total ice concentration | | 145 145 !! vt_i | - | Total ice vol. per unit area | m | 146 146 !! vt_s | - | Total snow vol. per unit ar. | m | … … 176 176 REAL(wp), PUBLIC :: ecc !: eccentricity of the elliptical yield curve 177 177 REAL(wp), PUBLIC :: ahi0 !: sea-ice hor. eddy diffusivity coeff. (m2/s) 178 REAL(wp), PUBLIC :: telast !: timescale for elastic waves (s) !SB 179 REAL(wp), PUBLIC :: alphaevp !: coeficient of the internal stresses !SB 178 REAL(wp), PUBLIC :: telast !: timescale for elastic waves (s) 179 REAL(wp), PUBLIC :: relast !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 180 REAL(wp), PUBLIC :: alphaevp !: coeficient of the internal stresses 180 181 REAL(wp), PUBLIC :: unit_fac = 1.e+09_wp !: conversion factor for ice / snow enthalpy 181 182 REAL(wp), PUBLIC :: hminrhg !: ice volume (a*h, in m) below which ice velocity is set to ocean velocity … … 247 248 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fhld !: heat flux from the lead used for bottom melting 248 249 249 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_snw !: Variation of snow mass over 1 time step [Kg/m2]250 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_ ice !: Variation of ice mass over 1 time step [Kg/m2]251 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_sub !: Variation of snow mass over 1 time step due to sublimation [Kg/m2]252 253 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_ sni !: snow ice growth254 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_ opw !: lateral ice growth255 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_ bog !: bottom ice growth256 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_ dyn !: dynamical ice growth257 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_ bom !: vertical bottom melt258 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_ sum !: vertical surface melt259 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_ res !: production (growth+melt) due to limupdate260 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_ spr !: snow precipitation on ice250 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_snw !: snow-ocean mass exchange over 1 time step [kg/m2] 251 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_spr !: snow precipitation on ice over 1 time step [kg/m2] 252 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_sub !: snow sublimation over 1 time step [kg/m2] 253 254 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_ice !: ice-ocean mass exchange over 1 time step [kg/m2] 255 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_sni !: snow ice growth component of wfx_ice [kg/m2] 256 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_opw !: lateral ice growth component of wfx_ice [kg/m2] 257 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_bog !: bottom ice growth component of wfx_ice [kg/m2] 258 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_dyn !: dynamical ice growth component of wfx_ice [kg/m2] 259 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_bom !: bottom melt component of wfx_ice [kg/m2] 260 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_sum !: surface melt component of wfx_ice [kg/m2] 261 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_res !: residual component of wfx_ice [kg/m2] 261 262 262 263 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx_bog !: salt flux due to ice growth/melt [PSU/m2/s] … … 323 324 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: smt_i !: mean sea ice salinity averaged over all categories [PSU] 324 325 325 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: at_i_typ !: total area contained in each ice type [m^2]326 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: vt_i_typ !: total volume contained in each ice type [m^3]327 328 326 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: t_s !: Snow temperatures [K] 329 327 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_s !: Snow ... 330 331 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_i_cat !: ! go to trash332 328 333 329 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: t_i !: ice temperatures [K] … … 350 346 !! * Old values of global variables 351 347 !!-------------------------------------------------------------------------- 352 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: old_v_s, old_v_i!: snow and ice volumes353 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: old_a_i, old_smv_i, old_oa_i !: ???354 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: old_e_s!: snow heat content355 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: old_e_i!: ice temperatures356 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: old_u_ice, old_v_ice !: ice velocity (gv6 and gv7)348 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_s_b, v_i_b !: snow and ice volumes 349 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_b, smv_i_b, oa_i_b !: 350 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_s_b !: snow heat content 351 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_i_b !: ice temperatures 352 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_ice_b, v_ice_b !: ice velocity 357 353 358 354 … … 377 373 !! * Ice thickness distribution variables 378 374 !!-------------------------------------------------------------------------- 379 ! REMOVE380 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ice_types !: Vector connecting types and categories381 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ice_cat_bounds !: Matrix containing the integer upper and382 ! ! lower boundaries of ice thickness categories383 ! REMOVE384 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ice_ncat_types !: nb of thickness categories in each ice type385 375 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hi_max !: Boundary of ice thickness categories in thickness space 386 376 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hi_mean !: Mean ice thickness in catgories 387 ! REMOVE388 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hi_max_typ !: Boundary of ice thickness categories in thickness space389 377 390 378 !!-------------------------------------------------------------------------- … … 475 463 & bv_i (jpi,jpj) , smt_i(jpi,jpj) , STAT=ierr(ii) ) 476 464 ii = ii + 1 477 ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , at_i_typ(jpi,jpj,jpm) ,&478 & e_s(jpi,jpj,nlay_s,jpl) , vt_i_typ(jpi,jpj,jpm) , e_i_cat(jpi,jpj,jpl) ,STAT=ierr(ii) )479 ii = ii + 1 480 ALLOCATE( t_i(jpi,jpj, jkmax,jpl) , e_i(jpi,jpj,jkmax,jpl) , s_i(jpi,jpj,jkmax,jpl) , STAT=ierr(ii) )465 ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , & 466 & e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) ) 467 ii = ii + 1 468 ALLOCATE( t_i(jpi,jpj,nlay_i+1,jpl) , e_i(jpi,jpj,nlay_i+1,jpl) , s_i(jpi,jpj,nlay_i+1,jpl) , STAT=ierr(ii) ) 481 469 482 470 ! * Moments for advection … … 494 482 & STAT=ierr(ii) ) 495 483 ii = ii + 1 496 ALLOCATE( sxe (jpi,jpj, jkmax,jpl) , sye (jpi,jpj,jkmax,jpl) , sxxe(jpi,jpj,jkmax,jpl) , &497 & syye(jpi,jpj, jkmax,jpl) , sxye(jpi,jpj,jkmax,jpl) , STAT=ierr(ii) )484 ALLOCATE( sxe (jpi,jpj,nlay_i+1,jpl) , sye (jpi,jpj,nlay_i+1,jpl) , sxxe(jpi,jpj,nlay_i+1,jpl) , & 485 & syye(jpi,jpj,nlay_i+1,jpl) , sxye(jpi,jpj,nlay_i+1,jpl) , STAT=ierr(ii) ) 498 486 499 487 ! * Old values of global variables 500 488 ii = ii + 1 501 ALLOCATE( old_v_s (jpi,jpj,jpl) , old_v_i (jpi,jpj,jpl) , old_e_s(jpi,jpj,nlay_s,jpl) , &502 & old_a_i (jpi,jpj,jpl) , old_smv_i(jpi,jpj,jpl) , old_e_i(jpi,jpj,jkmax,jpl) , &503 & o ld_oa_i(jpi,jpj,jpl) , &504 & old_u_ice(jpi,jpj) , old_v_ice(jpi,jpj) , STAT=ierr(ii) )489 ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) , & 490 & a_i_b (jpi,jpj,jpl) , smv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i+1 ,jpl) , & 491 & oa_i_b (jpi,jpj,jpl) , & 492 & u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , STAT=ierr(ii) ) 505 493 506 494 ! * Increment of global variables … … 512 500 & STAT=ierr(ii) ) 513 501 ii = ii + 1 514 ALLOCATE( d_e_s_thd(jpi,jpj,nlay_s,jpl) , d_e_i_thd(jpi,jpj, jkmax,jpl) , d_u_ice_dyn(jpi,jpj) , &515 & d_e_s_trp(jpi,jpj,nlay_s,jpl) , d_e_i_trp(jpi,jpj, jkmax,jpl) , d_v_ice_dyn(jpi,jpj) , STAT=ierr(ii) )502 ALLOCATE( d_e_s_thd(jpi,jpj,nlay_s,jpl) , d_e_i_thd(jpi,jpj,nlay_i+1,jpl) , d_u_ice_dyn(jpi,jpj) , & 503 & d_e_s_trp(jpi,jpj,nlay_s,jpl) , d_e_i_trp(jpi,jpj,nlay_i+1,jpl) , d_v_ice_dyn(jpi,jpj) , STAT=ierr(ii) ) 516 504 517 505 ! * Ice thickness distribution variables 518 506 ii = ii + 1 519 ALLOCATE( ice_types(jpl) , ice_cat_bounds(jpm,2) , ice_ncat_types (jpm) , & 520 & hi_max (0:jpl) , hi_mean(jpl) , hi_max_typ(0:jpl,jpm) , STAT=ierr(ii) ) 507 ALLOCATE( hi_max(0:jpl), hi_mean(jpl), STAT=ierr(ii) ) 521 508 522 509 ! * Ice diagnostics 523 510 ii = ii + 1 524 ALLOCATE( dv_dt_thd(jpi,jpj,jpl) ,&525 & izero (jpi,jpj,jpl) , diag_trp_vi(jpi,jpj) , diag_trp_vs(jpi,jpj), diag_trp_ei(jpi,jpj), diag_trp_es(jpi,jpj),&526 & diag_ heat_dhc(jpi,jpj), STAT=ierr(ii) )511 ALLOCATE( dv_dt_thd(jpi,jpj,jpl), izero (jpi,jpj,jpl), & 512 & diag_trp_vi(jpi,jpj), diag_trp_vs (jpi,jpj), diag_trp_ei(jpi,jpj), & 513 & diag_trp_es(jpi,jpj), diag_heat_dhc(jpi,jpj), STAT=ierr(ii) ) 527 514 528 515 ice_alloc = MAXVAL( ierr(:) ) -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
r4901 r4902 66 66 ! 67 67 ! ! adequation jpk versus ice/snow layers/categories 68 IF( jpl > jpk .OR. jpm > jpk .OR.&69 jkmax > jpk .OR. nlay_s > jpk ) CALL ctl_stop( 'STOP',&68 IF( jpl > jpk .OR. (nlay_i+1) > jpk .OR. nlay_s > jpk ) & 69 & CALL ctl_stop( 'STOP', & 70 70 & 'ice_init: the 3rd dimension of workspace arrays is too small.', & 71 71 & 'use more ocean levels or less ice/snow layers/categories.' ) … … 174 174 !! limistate (only) and is changed to 99 m in ice_init 175 175 !!------------------------------------------------------------------ 176 INTEGER :: jl , jm! dummy loop index176 INTEGER :: jl ! dummy loop index 177 177 REAL(wp) :: zc1, zc2, zc3, zx1 ! local scalars 178 178 !!------------------------------------------------------------------ … … 185 185 ! 1) Ice thickness distribution parameters initialization 186 186 !------------------------------------------------------------------------------! 187 188 !- Types boundaries (integer)189 !----------------------------190 ice_cat_bounds(1,1) = 1191 ice_cat_bounds(1,2) = jpl192 193 !- Number of ice thickness categories in each ice type194 DO jm = 1, jpm195 ice_ncat_types(jm) = ice_cat_bounds(jm,2) - ice_cat_bounds(jm,1) + 1196 END DO197 198 !- Make the correspondence between thickness categories and ice types199 !---------------------------------------------------------------------200 DO jm = 1, jpm !over types201 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) !over thickness categories202 ice_types(jl) = jm203 END DO204 END DO205 206 187 IF(lwp) THEN 207 WRITE(numout,*) ' Number of ice types jpm = ', jpm208 188 WRITE(numout,*) ' Number of ice categories jpl = ', jpl 209 DO jm = 1, jpm210 WRITE(numout,*) ' Ice type ', jm211 WRITE(numout,*) ' Number of thickness categories ', ice_ncat_types(jm)212 WRITE(numout,*) ' Thickness category boundaries ', ice_cat_bounds(jm,1:2)213 END DO214 WRITE(numout,*) 'Ice type vector', ice_types(1:jpl)215 WRITE(numout,*)216 189 ENDIF 217 190 … … 219 192 !---------------------------------- 220 193 hi_max(:) = 0._wp 221 hi_max_typ(:,:) = 0._wp 222 223 !- Type 1 - undeformed ice 224 zc1 = 3._wp / REAL( ice_cat_bounds(1,2) - ice_cat_bounds(1,1) + 1 , wp ) 194 195 zc1 = 3._wp / REAL( jpl, wp ) 225 196 zc2 = 10._wp * zc1 226 197 zc3 = 3._wp 227 198 228 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)229 zx1 = REAL( jl-1 , wp ) / REAL( ice_cat_bounds(1,2) - ice_cat_bounds(1,1) + 1, wp )199 DO jl = 1, jpl 200 zx1 = REAL( jl-1, wp ) / REAL( jpl, wp ) 230 201 hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) ) 231 202 END DO 232 203 233 !- Fill in the hi_max_typ vector, useful in other circumstances 234 ! Tricky trick: hi_max_typ is actually not used in the code and will be removed in a 235 ! next flyspray at this time, the tricky trick will also be removed (Martin, march 08) 236 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 237 hi_max_typ(jl,1) = hi_max(jl) 238 END DO 239 240 IF(lwp) WRITE(numout,*) ' Thickness category boundaries independently of ice type ' 204 IF(lwp) WRITE(numout,*) ' Thickness category boundaries ' 241 205 IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl) 242 206 243 IF(lwp) WRITE(numout,*) ' Thickness category boundaries inside ice types '244 IF(lwp) THEN245 DO jm = 1, jpm246 WRITE(numout,*) ' Type number ', jm247 WRITE(numout,*) ' hi_max_typ : ', hi_max_typ(0:ice_ncat_types(jm),jm)248 END DO249 ENDIF250 207 ! 251 208 DO jl = 1, jpl -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
r4900 r4902 73 73 !! ** Method : Arithmetics 74 74 !!--------------------------------------------------------------------- 75 INTEGER , INTENT(in ) :: ksum !: number of categories76 INTEGER , INTENT(in ) :: klay !: number of vertical layers77 REAL(wp), DIMENSION(jpi,jpj, jkmax,jpl), INTENT(in ) :: pin !: input field78 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pout !: output field75 INTEGER , INTENT(in ) :: ksum !: number of categories 76 INTEGER , INTENT(in ) :: klay !: number of vertical layers 77 REAL(wp), DIMENSION(jpi,jpj,nlay_i+1,jpl), INTENT(in ) :: pin !: input field 78 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pout !: output field 79 79 ! 80 80 INTEGER :: jk, jl ! dummy loop indices … … 175 175 zei_b = glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 176 176 zfw_b = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 177 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) ) * area(:,:) * tms(:,:) ) 177 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) & 178 & ) * area(:,:) * tms(:,:) ) 178 179 zfs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 179 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 180 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) & 181 & ) * area(:,:) * tms(:,:) ) 180 182 zft_b = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 181 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 183 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 184 & ) * area(:,:) / unit_fac * tms(:,:) ) 182 185 183 186 ELSEIF( icount == 1 ) THEN 184 187 185 188 zfs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 186 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zfs_b 189 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) & 190 & ) * area(:,:) * tms(:,:) ) - zfs_b 187 191 zfw = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 188 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) ) * area(:,:) * tms(:,:) ) - zfw_b 192 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) & 193 & ) * area(:,:) * tms(:,:) ) - zfw_b 189 194 zft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 190 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zft_b 195 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 196 & ) * area(:,:) / unit_fac * tms(:,:) ) - zft_b 191 197 192 198 zvi = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zvi_b ) * r1_rdtice - zfw -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
r4900 r4902 59 59 !! 60 60 REAL(dp) :: zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc 61 REAL(dp) :: zbg_sfx, zbg_sfx_bri, zbg_sfx_bog, zbg_sfx_bom, zbg_sfx_sum, zbg_sfx_sni, zbg_sfx_opw, zbg_sfx_res, zbg_sfx_dyn 61 REAL(dp) :: zbg_sfx, zbg_sfx_bri, zbg_sfx_bog, zbg_sfx_bom, zbg_sfx_sum, zbg_sfx_sni, & 62 & zbg_sfx_opw, zbg_sfx_res, zbg_sfx_dyn 62 63 REAL(dp) :: zbg_vfx, zbg_vfx_bog, zbg_vfx_opw, zbg_vfx_sni, zbg_vfx_dyn 63 64 REAL(dp) :: zbg_vfx_bom, zbg_vfx_sum, zbg_vfx_res, zbg_vfx_spr, zbg_vfx_snw, zbg_vfx_sub -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r4900 r4902 83 83 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 84 84 85 old_u_ice(:,:) = u_ice(:,:) * tmu(:,:)86 old_v_ice(:,:) = v_ice(:,:) * tmv(:,:)85 u_ice_b(:,:) = u_ice(:,:) * tmu(:,:) 86 v_ice_b(:,:) = v_ice(:,:) * tmv(:,:) 87 87 88 88 ! Rheology (ice dynamics) … … 243 243 NAMELIST/namicedyn/ epsd, om, cw, angvg, pstar, & 244 244 & c_rhg, creepl, ecc, ahi0, & 245 & nevp, telast, alphaevp, hminrhg245 & nevp, relast, alphaevp, hminrhg 246 246 !!------------------------------------------------------------------- 247 247 … … 269 269 WRITE(numout,*) ' horizontal diffusivity coeff. for sea-ice ahi0 = ', ahi0 270 270 WRITE(numout,*) ' number of iterations for subcycling nevp = ', nevp 271 WRITE(numout,*) ' timescale for elastic waves telast = ', telast271 WRITE(numout,*) ' ratio of elastic timescale over ice time step relast = ', relast 272 272 WRITE(numout,*) ' coefficient for the solution of int. stresses alphaevp = ', alphaevp 273 273 WRITE(numout,*) ' min ice thickness for rheology calculations hminrhg = ', hminrhg … … 287 287 pstarh = pstar * 0.5_wp 288 288 289 ! elastic damping 290 telast = relast * rdt_ice 291 289 292 ! Diffusion coefficients. 290 293 ahiu(:,:) = ahi0 * umask(:,:,1) -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r4901 r4902 128 128 DO jj = 1, jpj ! ice if sst <= t-freez + ttest 129 129 DO ji = 1, jpi 130 IF( ( tsn(ji,jj,1,jp_tem) - ( t_bo(ji,jj) - rt0 ) ) * tms(ji,jj) >= thres_sst ) THEN ; zswitch(ji,jj) = 0._wp * tms(ji,jj) ! no ice 131 ELSE ; zswitch(ji,jj) = 1._wp * tms(ji,jj) ! ice 130 IF( ( tsn(ji,jj,1,jp_tem) - ( t_bo(ji,jj) - rt0 ) ) * tms(ji,jj) >= thres_sst ) THEN 131 zswitch(ji,jj) = 0._wp * tms(ji,jj) ! no ice 132 ELSE 133 zswitch(ji,jj) = 1._wp * tms(ji,jj) ! ice 132 134 ENDIF 133 135 END DO -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r4901 r4902 692 692 693 693 IF( partfun_swi == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975) 694 DO jl = 0, ice_cat_bounds(1,2) ! only undeformed ice participates694 DO jl = 0, jpl 695 695 DO jj = 1, jpj 696 696 DO ji = 1, jpi … … 715 715 Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 716 716 END DO !jl 717 DO jl = 0, ice_cat_bounds(1,2)717 DO jl = 0, jpl 718 718 athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 719 719 END DO … … 897 897 CALL wrk_alloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 898 898 CALL wrk_alloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 899 CALL wrk_alloc( jpi, jpj, jkmax, eirft, erdg1, erdg2, ersw )900 CALL wrk_alloc( jpi, jpj, jkmax, jpl, eicen_init )899 CALL wrk_alloc( jpi, jpj, nlay_i+1, eirft, erdg1, erdg2, ersw ) 900 CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, eicen_init ) 901 901 902 902 ! Conservation check … … 1193 1193 !------------------------------------------------------------------------------- 1194 1194 ! jl1 looping 1-jpl 1195 DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2)1195 DO jl2 = 1, jpl 1196 1196 ! over categories to which ridged ice is transferred 1197 1197 !CDIR NODEP … … 1238 1238 END DO ! jl2 (new ridges) 1239 1239 1240 DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2)1240 DO jl2 = 1, jpl 1241 1241 1242 1242 !CDIR NODEP … … 1302 1302 CALL wrk_dealloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 1303 1303 CALL wrk_dealloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 1304 CALL wrk_dealloc( jpi, jpj, jkmax, eirft, erdg1, erdg2, ersw )1305 CALL wrk_dealloc( jpi, jpj, jkmax, jpl, eicen_init )1304 CALL wrk_dealloc( jpi, jpj, nlay_i+1, eirft, erdg1, erdg2, ersw ) 1305 CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, eicen_init ) 1306 1306 ! 1307 1307 END SUBROUTINE lim_itd_me_ridgeshift -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r4900 r4902 6 6 !! History : - ! (W. H. Lipscomb and E.C. Hunke) CICE (c) original code 7 7 !! 3.0 ! 2005-12 (M. Vancoppenolle) adaptation to LIM-3 8 !! - ! 2006-06 (M. Vancoppenolle) adaptation to include salt, age and types8 !! - ! 2006-06 (M. Vancoppenolle) adaptation to include salt, age 9 9 !! - ! 2007-04 (M. Vancoppenolle) Mass conservation checked 10 10 !!---------------------------------------------------------------------- … … 66 66 INTEGER, INTENT(in) :: kt ! time step index 67 67 ! 68 INTEGER :: ji, jj, jk, jl, ja, jm, jbnd1, jbnd2 ! ice typesdummy loop index68 INTEGER :: ji, jj, jk, jl ! dummy loop index 69 69 ! 70 70 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b … … 86 86 ! Given thermodynamic growth rates, transport ice between 87 87 ! thickness categories. 88 DO jm = 1, jpm 89 jbnd1 = ice_cat_bounds(jm,1) 90 jbnd2 = ice_cat_bounds(jm,2) 91 IF( ice_ncat_types(jm) > 1 ) CALL lim_itd_th_rem( jbnd1, jbnd2, jm, kt ) 92 END DO 88 IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt ) 93 89 ! 94 90 CALL lim_var_glo2eqv ! only for info … … 123 119 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_th : sm_i : ') 124 120 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_th : smv_i : ') 125 DO j a= 1, nlay_i121 DO jk = 1, nlay_i 126 122 CALL prt_ctl_info(' ') 127 CALL prt_ctl_info(' - Layer : ', ivar1=j a)123 CALL prt_ctl_info(' - Layer : ', ivar1=jk) 128 124 CALL prt_ctl_info(' ~~~~~~~') 129 CALL prt_ctl(tab2d_1=t_i(:,:,j a,jl) , clinfo1= ' lim_itd_th : t_i : ')130 CALL prt_ctl(tab2d_1=e_i(:,:,j a,jl) , clinfo1= ' lim_itd_th : e_i : ')125 CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : t_i : ') 126 CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : e_i : ') 131 127 END DO 132 128 END DO … … 140 136 ! 141 137 142 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, ntyp,kt )138 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt ) 143 139 !!------------------------------------------------------------------ 144 140 !! *** ROUTINE lim_itd_th_rem *** … … 153 149 INTEGER , INTENT (in) :: klbnd ! Start thickness category index point 154 150 INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied 155 INTEGER , INTENT (in) :: ntyp ! Number of the type used156 151 INTEGER , INTENT (in) :: kt ! Ocean time step 157 152 ! … … 171 166 REAL(wp), POINTER, DIMENSION(:,:,:) :: hL ! left boundary for the ITD for each thickness 172 167 REAL(wp), POINTER, DIMENSION(:,:,:) :: hR ! left boundary for the ITD for each thickness 173 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_ o! old ice thickness168 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_b ! old ice thickness 174 169 REAL(wp), POINTER, DIMENSION(:,:,:) :: dummy_es 175 170 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! local increment of ice area and volume … … 189 184 CALL wrk_alloc( jpi,jpj, zremap_flag ) ! integer 190 185 CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) ! integer 191 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_ o, dummy_es )186 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 192 187 CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice ) 193 188 CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) … … 218 213 WRITE(numout,*) ' klbnd : ', klbnd 219 214 WRITE(numout,*) ' kubnd : ', kubnd 220 WRITE(numout,*) ' ntyp : ', ntyp221 215 ENDIF 222 216 … … 227 221 zindb = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 228 222 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb 229 zindb = 1.0 - MAX( 0.0, SIGN( 1.0, - old_a_i(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes230 zht_i_ o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX( old_a_i(ji,jj,jl), epsi10 ) * zindb231 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_ o(ji,jj,jl)223 zindb = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 224 zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * zindb 225 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) 232 226 END DO 233 227 END DO … … 274 268 ! 275 269 zhbnew(ii,ij,jl) = hi_max(jl) 276 IF ( old_a_i(ii,ij,jl) > epsi10 .AND. old_a_i(ii,ij,jl+1) > epsi10 ) THEN270 IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 277 271 !interpolate between adjacent category growth rates 278 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_ o(ii,ij,jl+1) - zht_i_o(ii,ij,jl) )279 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_ o(ii,ij,jl) )280 ELSEIF ( old_a_i(ii,ij,jl) > epsi10) THEN272 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 273 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 274 ELSEIF ( a_i_b(ii,ij,jl) > epsi10) THEN 281 275 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 282 ELSEIF ( old_a_i(ii,ij,jl+1) > epsi10) THEN276 ELSEIF ( a_i_b(ii,ij,jl+1) > epsi10) THEN 283 277 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 284 278 ENDIF … … 321 315 DO jj = 1, jpj 322 316 DO ji = 1, jpi 323 zhb0(ji,jj) = hi_max _typ(0,ntyp) ! 0eme324 zhb1(ji,jj) = hi_max _typ(1,ntyp) ! 1er317 zhb0(ji,jj) = hi_max(0) ! 0eme 318 zhb1(ji,jj) = hi_max(1) ! 1er 325 319 326 320 zhbnew(ji,jj,klbnd-1) = 0._wp … … 343 337 !----------------------------------------------------------------------------------------------- 344 338 !- 7.1 g(h) for category 1 at start of time step 345 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_ o(:,:,klbnd), &339 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), & 346 340 & g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 347 341 & hR(:,:,klbnd), zremap_flag ) … … 368 362 ! Constrain new thickness <= ht_i 369 363 zdamax = a_i(ii,ij,klbnd) * & 370 (1.0 - ht_i(ii,ij,klbnd)/zht_i_ o(ii,ij,klbnd)) ! zdamax > 0364 (1.0 - ht_i(ii,ij,klbnd)/zht_i_b(ii,ij,klbnd)) ! zdamax > 0 371 365 !ice area lost due to melting of thin ice 372 366 zda0 = MIN(zda0, zdamax) … … 382 376 ELSE ! if ice accretion 383 377 ! ji, a_i > epsi10; zdh0 > 0 384 IF ( ntyp .EQ. 1 )zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd))378 zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 385 379 ! zhbnew was 0, and is shifted to the right to account for thin ice 386 380 ! growth in openwater (F0 = f1) 387 IF ( ntyp .NE. 1 ) zhbnew(ii,ij,0) = 0388 ! in other types there is389 ! no open water growth (F0 = 0)390 381 ENDIF ! zdh0 391 382 … … 493 484 CALL wrk_dealloc( jpi,jpj, zremap_flag ) ! integer 494 485 CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) ! integer 495 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_ o, dummy_es )486 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 496 487 CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice ) 497 488 CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) … … 839 830 840 831 841 SUBROUTINE lim_itd_th_reb( klbnd, kubnd , ntyp)832 SUBROUTINE lim_itd_th_reb( klbnd, kubnd ) 842 833 !!------------------------------------------------------------------ 843 834 !! *** ROUTINE lim_itd_th_reb *** … … 849 840 INTEGER , INTENT (in) :: klbnd ! Start thickness category index point 850 841 INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied 851 INTEGER , INTENT (in) :: ntyp ! number of the ice type involved in the rebinning process852 842 ! 853 843 INTEGER :: ji,jj, jl ! dummy loop indices … … 889 879 890 880 !------------------------------------------------------------------------------ 891 ! 2) Make sure thickness of cat klbnd is at least hi_max _typ(klbnd)881 ! 2) Make sure thickness of cat klbnd is at least hi_max(klbnd) 892 882 !------------------------------------------------------------------------------ 893 883 DO jj = 1, jpj 894 884 DO ji = 1, jpi 895 885 IF( a_i(ji,jj,klbnd) > epsi10 ) THEN 896 IF( ht_i(ji,jj,klbnd) <= hi_max _typ(0,ntyp) .AND. hi_max_typ(0,ntyp) > 0._wp ) THEN897 a_i(ji,jj,klbnd) = v_i(ji,jj,klbnd) / hi_max _typ(0,ntyp)898 ht_i(ji,jj,klbnd) = hi_max _typ(0,ntyp)886 IF( ht_i(ji,jj,klbnd) <= hi_max(0) .AND. hi_max(0) > 0._wp ) THEN 887 a_i(ji,jj,klbnd) = v_i(ji,jj,klbnd) / hi_max(0) 888 ht_i(ji,jj,klbnd) = hi_max(0) 899 889 ENDIF 900 890 ENDIF -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r4901 r4902 118 118 CALL iom_put( "qsr_oce" , qsr(:,:) * pfrld(:,:) ) ! solar flux at ocean surface 119 119 CALL iom_put( "qns_oce" , qns(:,:) * pfrld(:,:) ) ! non-solar flux at ocean surface 120 CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) ) ! solar flux at ice surface121 CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) ) ! non-solar flux at ice surface122 CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) ) ! solar flux transmitted thru ice120 CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux at ice surface 121 CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! non-solar flux at ice surface 122 CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux transmitted thru ice 123 123 CALL iom_put( "qt_oce" , ( qsr(:,:) + qns(:,:) ) * pfrld(:,:) ) 124 CALL iom_put( "qt_ice" , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) * old_a_i(:,:,:), dim=3 ) )124 CALL iom_put( "qt_ice" , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) * a_i_b(:,:,:), dim=3 ) ) 125 125 126 126 ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) … … 139 139 zfcm1 = qsr_tot(ji,jj) 140 140 DO jl = 1, jpl 141 zfcm1 = zfcm1 + ( ftr_ice(ji,jj,jl) - qsr_ice(ji,jj,jl) ) * old_a_i(ji,jj,jl)141 zfcm1 = zfcm1 + ( ftr_ice(ji,jj,jl) - qsr_ice(ji,jj,jl) ) * a_i_b(ji,jj,jl) 142 142 END DO 143 143 ELSE … … 145 145 zfcm1 = pfrld(ji,jj) * qsr(ji,jj) 146 146 DO jl = 1, jpl 147 zfcm1 = zfcm1 + old_a_i(ji,jj,jl) * ftr_ice(ji,jj,jl)147 zfcm1 = zfcm1 + a_i_b(ji,jj,jl) * ftr_ice(ji,jj,jl) 148 148 END DO 149 149 ENDIF … … 181 181 182 182 ! mass flux from ice/ocean 183 wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) 183 wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) & 184 + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) 184 185 185 186 ! mass flux at the ocean/ice interface 186 187 fmmflx(ji,jj) = - wfx_ice(ji,jj) * rdt_ice ! F/M mass flux save at least for biogeochemical model 187 emp(ji,jj) = zemp - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_sub(ji,jj)! mass flux + F/M mass flux (always ice/ocean mass exchange)188 emp(ji,jj) = zemp - wfx_ice(ji,jj) - wfx_snw(ji,jj) ! mass flux + F/M mass flux (always ice/ocean mass exchange) 188 189 189 190 END DO … … 193 194 ! salt flux at the ocean surface ! 194 195 !------------------------------------------! 195 sfx(:,:) = sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:) 196 sfx(:,:) = sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) & 197 & + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:) 196 198 197 199 !-------------------------------------------------------------! -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r4901 r4902 22 22 USE phycst ! physical constants 23 23 USE dom_oce ! ocean space and time domain variables 24 USE oce , ONLY : fraqsr_1lev 24 USE oce , ONLY : fraqsr_1lev 25 25 USE ice ! LIM: sea-ice variables 26 26 USE par_ice ! LIM: sea-ice parameters … … 153 153 DO jj = 1, jpj 154 154 DO ji = 1, jpi 155 zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * old_a_i(ji,jj,jl)156 zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * old_a_i(ji,jj,jl)155 zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * a_i_b(ji,jj,jl) 156 zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * a_i_b(ji,jj,jl) 157 157 END DO 158 158 END DO … … 282 282 !------------------------- 283 283 284 CALL tab_2d_1d( nbpb, at_i_ b(1:nbpb), at_i , jpi, jpj, npb(1:nbpb) )285 CALL tab_2d_1d( nbpb, a_i_ b(1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )286 CALL tab_2d_1d( nbpb, ht_i_ b(1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )287 CALL tab_2d_1d( nbpb, ht_s_ b(1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) )288 289 CALL tab_2d_1d( nbpb, t_su_ b(1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) )290 CALL tab_2d_1d( nbpb, sm_i_ b(1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )284 CALL tab_2d_1d( nbpb, at_i_1d (1:nbpb), at_i , jpi, jpj, npb(1:nbpb) ) 285 CALL tab_2d_1d( nbpb, a_i_1d (1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 286 CALL tab_2d_1d( nbpb, ht_i_1d (1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 287 CALL tab_2d_1d( nbpb, ht_s_1d (1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 288 289 CALL tab_2d_1d( nbpb, t_su_1d (1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 290 CALL tab_2d_1d( nbpb, sm_i_1d (1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 291 291 DO jk = 1, nlay_s 292 CALL tab_2d_1d( nbpb, t_s_ b(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )293 CALL tab_2d_1d( nbpb, q_s_ b(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )292 CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 293 CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 294 294 END DO 295 295 DO jk = 1, nlay_i 296 CALL tab_2d_1d( nbpb, t_i_ b(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )297 CALL tab_2d_1d( nbpb, q_i_ b(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )298 CALL tab_2d_1d( nbpb, s_i_ b(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )296 CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 297 CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 298 CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 299 299 END DO 300 300 … … 310 310 ENDIF 311 311 CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 312 CALL tab_2d_1d( nbpb, t_bo_ b(1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) )312 CALL tab_2d_1d( nbpb, t_bo_1d (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) ) 313 313 CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip , jpi, jpj, npb(1:nbpb) ) 314 314 CALL tab_2d_1d( nbpb, fhtur_1d (1:nbpb), fhtur , jpi, jpj, npb(1:nbpb) ) … … 361 361 362 362 ! --- Ice enthalpy remapping --- ! 363 CALL lim_thd_ent( 1, nbpb, q_i_ b(1:nbpb,:) )363 CALL lim_thd_ent( 1, nbpb, q_i_1d(1:nbpb,:) ) 364 364 365 365 !---------------------------------! … … 377 377 !-------------------------------- 378 378 379 CALL tab_1d_2d( nbpb, at_i , npb, at_i_ b(1:nbpb) , jpi, jpj )380 CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_ b(1:nbpb) , jpi, jpj )381 CALL tab_1d_2d( nbpb, ht_s(:,:,jl) , npb, ht_s_ b(1:nbpb) , jpi, jpj )382 CALL tab_1d_2d( nbpb, a_i (:,:,jl) , npb, a_i_ b(1:nbpb) , jpi, jpj )383 CALL tab_1d_2d( nbpb, t_su(:,:,jl) , npb, t_su_ b(1:nbpb) , jpi, jpj )384 CALL tab_1d_2d( nbpb, sm_i(:,:,jl) , npb, sm_i_ b(1:nbpb) , jpi, jpj )379 CALL tab_1d_2d( nbpb, at_i , npb, at_i_1d (1:nbpb) , jpi, jpj ) 380 CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_1d (1:nbpb) , jpi, jpj ) 381 CALL tab_1d_2d( nbpb, ht_s(:,:,jl) , npb, ht_s_1d (1:nbpb) , jpi, jpj ) 382 CALL tab_1d_2d( nbpb, a_i (:,:,jl) , npb, a_i_1d (1:nbpb) , jpi, jpj ) 383 CALL tab_1d_2d( nbpb, t_su(:,:,jl) , npb, t_su_1d (1:nbpb) , jpi, jpj ) 384 CALL tab_1d_2d( nbpb, sm_i(:,:,jl) , npb, sm_i_1d (1:nbpb) , jpi, jpj ) 385 385 DO jk = 1, nlay_s 386 CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_ b(1:nbpb,jk), jpi, jpj)387 CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_ b(1:nbpb,jk), jpi, jpj)386 CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d (1:nbpb,jk), jpi, jpj) 387 CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d (1:nbpb,jk), jpi, jpj) 388 388 END DO 389 389 DO jk = 1, nlay_i 390 CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_ b(1:nbpb,jk), jpi, jpj)391 CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_ b(1:nbpb,jk), jpi, jpj)392 CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_ b(1:nbpb,jk), jpi, jpj)390 CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d (1:nbpb,jk), jpi, jpj) 391 CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d (1:nbpb,jk), jpi, jpj) 392 CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d (1:nbpb,jk), jpi, jpj) 393 393 END DO 394 394 CALL tab_1d_2d( nbpb, qlead , npb, qlead_1d (1:nbpb) , jpi, jpj ) … … 532 532 DO jk = 1, nlay_i 533 533 DO ji = kideb, kiut 534 ztmelts = -tmut * s_i_ b(ji,jk) + rtt534 ztmelts = -tmut * s_i_1d(ji,jk) + rtt 535 535 ! Conversion q(S,T) -> T (second order equation) 536 536 zaaa = cpic 537 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_ b(ji,jk) / rhoic - lfus537 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_1d(ji,jk) / rhoic - lfus 538 538 zccc = lfus * ( ztmelts - rtt ) 539 539 zdiscrim = SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) 540 t_i_ b(ji,jk) = rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa )540 t_i_1d(ji,jk) = rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) 541 541 542 542 ! mask temperature 543 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_ b(ji) ) )544 t_i_ b(ji,jk) = zswitch * t_i_b(ji,jk) + ( 1._wp - zswitch ) * rtt543 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 544 t_i_1d(ji,jk) = zswitch * t_i_1d(ji,jk) + ( 1._wp - zswitch ) * rtt 545 545 END DO 546 546 END DO -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r4901 r4902 129 129 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 130 130 CALL wrk_alloc( jpij, zintermelt ) 131 CALL wrk_alloc( jpij, jkmax, zdeltah, zh_i )131 CALL wrk_alloc( jpij, nlay_i+1, zdeltah, zh_i ) 132 132 CALL wrk_alloc( jpij, icount ) 133 133 … … 155 155 DO jk = 1, nlay_i 156 156 DO ji = kideb, kiut 157 h_i_old (ji,jk) = ht_i_ b(ji) / REAL( nlay_i )158 qh_i_old(ji,jk) = q_i_ b(ji,jk) * h_i_old(ji,jk)157 h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 158 qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk) 159 159 ENDDO 160 160 ENDDO … … 165 165 ! 166 166 DO ji = kideb, kiut 167 zinda = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )168 ztmelts = zinda * rtt + ( 1._wp - zinda ) * rtt167 zinda = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 168 ztmelts = zinda * rtt + ( 1._wp - zinda ) * rtt 169 169 170 170 zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 171 171 zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 172 172 173 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_ b(ji) - ztmelts ) )173 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts ) ) 174 174 zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 175 175 END DO … … 181 181 !------------------------------------------------------------------------------! 182 182 DO ji = kideb, kiut 183 IF( t_s_ b(ji,1) > rtt ) THEN !!! Internal melting183 IF( t_s_1d(ji,1) > rtt ) THEN !!! Internal melting 184 184 ! Contribution to heat flux to the ocean [W.m-2], < 0 185 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_ b(ji,1) * ht_s_b(ji) * a_i_b(ji) * r1_rdtice185 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 186 186 ! Contribution to mass flux 187 wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_ b(ji) * a_i_b(ji) * r1_rdtice187 wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 188 188 ! updates 189 ht_s_ b(ji) = 0._wp190 q_s_ b(ji,1) = 0._wp191 t_s_ b(ji,1) = rtt189 ht_s_1d(ji) = 0._wp 190 q_s_1d (ji,1) = 0._wp 191 t_s_1d (ji,1) = rtt 192 192 END IF 193 193 END DO … … 198 198 ! 199 199 DO ji = kideb, kiut 200 zh_s(ji) = ht_s_ b(ji) / REAL( nlay_s )200 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 201 201 END DO 202 202 ! 203 203 DO jk = 1, nlay_s 204 204 DO ji = kideb, kiut 205 zqh_s(ji) = zqh_s(ji) + q_s_ b(ji,jk) * zh_s(ji)205 zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji) 206 206 END DO 207 207 END DO … … 209 209 DO jk = 1, nlay_i 210 210 DO ji = kideb, kiut 211 zh_i(ji,jk) = ht_i_ b(ji) / REAL( nlay_i )212 zqh_i(ji) = zqh_i(ji) + q_i_ b(ji,jk) * zh_i(ji,jk)211 zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 212 zqh_i(ji) = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk) 213 213 END DO 214 214 END DO … … 237 237 !----------- 238 238 ! thickness change 239 zcoeff = ( 1._wp - ( 1._wp - at_i_ b(ji) )**betas ) / at_i_b(ji)239 zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**betas ) / at_i_1d(ji) 240 240 zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn 241 241 ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) … … 243 243 IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 244 244 ! heat flux from snow precip (>0, W.m-2) 245 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_ b(ji) * zqprec(ji) * r1_rdtice245 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 246 246 ! mass flux, <0 247 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_ b(ji) * zdh_s_pre(ji) * r1_rdtice247 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice 248 248 ! update thickness 249 ht_s_ b (ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_pre(ji) )249 ht_s_1d (ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) 250 250 251 251 !--------------------- … … 258 258 zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting 259 259 ! heat used to melt snow (W.m-2, >0) 260 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_ b(ji) * zqprec(ji) * r1_rdtice260 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 261 261 ! snow melting only = water into the ocean (then without snow precip), >0 262 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_ b(ji) * zdh_s_mel(ji) * r1_rdtice262 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_mel(ji) * r1_rdtice 263 263 264 264 ! updates available heat + thickness 265 265 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) ) 266 ht_s_ b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_mel(ji) )267 zh_s (ji) = ht_s_ b(ji) / REAL( nlay_s )266 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) ) 267 zh_s (ji) = ht_s_1d(ji) / REAL( nlay_s ) 268 268 269 269 ENDIF … … 275 275 DO ji = kideb, kiut 276 276 ! thickness change 277 zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_ b(ji) ) )278 zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_ b(ji,jk) + epsi20 ) )279 zdeltah (ji,jk) = - zindh * zindq * zq_su(ji) / MAX( q_s_ b(ji,jk), epsi20 )277 zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 278 zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_1d(ji,jk) + epsi20 ) ) 279 zdeltah (ji,jk) = - zindh * zindq * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 ) 280 280 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting 281 281 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 282 282 ! heat used to melt snow(W.m-2, >0) 283 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_ b(ji) * q_s_b(ji,jk) * r1_rdtice283 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice 284 284 ! snow melting only = water into the ocean (then without snow precip) 285 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice285 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 286 286 287 287 ! updates available heat + thickness 288 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_ b(ji,jk) )289 ht_s_ b(ji) = MAX( 0._wp , ht_s_b(ji) + zdeltah(ji,jk) )288 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 289 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 290 290 291 291 END DO … … 304 304 ! forced mode: snow thickness change due to sublimation 305 305 DO ji = kideb, kiut 306 zdh_s_sub(ji) = MAX( - ht_s_ b(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice )306 zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) 307 307 ! Heat flux by sublimation [W.m-2], < 0 308 308 ! sublimate first snow that had fallen, then pre-existing snow 309 309 zcoeff = ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) * zqprec(ji) + & 310 & ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_ b(ji,1) ) &311 & * a_i_ b(ji) * r1_rdtice310 & ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_1d(ji,1) ) & 311 & * a_i_1d(ji) * r1_rdtice 312 312 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff 313 313 ! Mass flux by sublimation 314 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_ b(ji) * zdh_s_sub(ji) * r1_rdtice314 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 315 315 ! new snow thickness 316 ht_s_ b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_sub(ji) )316 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 317 317 END DO 318 318 ENDIF … … 321 321 DO ji = kideb, kiut 322 322 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 323 zh_s(ji) = ht_s_ b(ji) / REAL( nlay_s )323 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 324 324 END DO ! ji 325 325 … … 331 331 DO jk = 1, nlay_s 332 332 DO ji = kideb,kiut 333 zindh = MAX( 0._wp , SIGN( 1._wp, - ht_s_ b(ji) + epsi20 ) )334 q_s_ b(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_b(ji), epsi20 ) * &333 zindh = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) + epsi20 ) ) 334 q_s_1d(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_1d(ji), epsi20 ) * & 335 335 & ( ( MAX( 0._wp, dh_s_tot(ji) ) ) * zqprec(ji) + & 336 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_ b(ji) ) * rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) )337 zq_s(ji) = zq_s(ji) + q_s_ b(ji,jk)336 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) ) 337 zq_s(ji) = zq_s(ji) + q_s_1d(ji,jk) 338 338 END DO 339 339 END DO … … 345 345 DO jk = 1, nlay_i 346 346 DO ji = kideb, kiut 347 zEi = - q_i_ b(ji,jk) / rhoic ! Specific enthalpy of layer k [J/kg, <0]348 349 ztmelts = - tmut * s_i_ b(ji,jk) + rtt ! Melting point of layer k [K]347 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of layer k [J/kg, <0] 348 349 ztmelts = - tmut * s_i_1d(ji,jk) + rtt ! Melting point of layer k [K] 350 350 351 351 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] … … 367 367 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 368 368 369 ! Contribution to salt flux (clem: using sm_i_ b and not s_i_b(jk) is ok)370 sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_ b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice369 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 370 sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 371 371 372 372 ! Contribution to heat flux [W.m-2], < 0 373 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice373 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 374 374 375 375 ! Total heat flux used in this process [W.m-2], > 0 376 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_ b(ji) * zdE * r1_rdtice376 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 377 377 378 378 ! Contribution to mass flux 379 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice379 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 380 380 381 381 ! record which layers have disappeared (for bottom melting) … … 387 387 388 388 ! update heat content (J.m-2) and layer thickness 389 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_ b(ji,jk)389 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 390 390 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 391 391 END DO … … 393 393 ! update ice thickness 394 394 DO ji = kideb, kiut 395 ht_i_ b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_surf(ji) )395 ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) ) 396 396 END DO 397 397 … … 423 423 !clem debug. Just to be sure that enthalpy at nlay_i+1 is null 424 424 DO ji = kideb, kiut 425 q_i_ b(ji,nlay_i+1) = 0._wp425 q_i_1d(ji,nlay_i+1) = 0._wp 426 426 END DO 427 427 … … 445 445 446 446 s_i_new(ji) = zswitch_sal * zfracs * sss_m(ii,ij) & ! New ice salinity 447 + ( 1. - zswitch_sal ) * sm_i_ b(ji)447 + ( 1. - zswitch_sal ) * sm_i_1d(ji) 448 448 ! New ice growth 449 449 ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K) 450 450 451 zt_i_new = zswitch_sal * t_bo_ b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)451 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 452 452 453 453 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) … … 455 455 & + rcp * ( ztmelts-rtt ) 456 456 457 zEw = rcp * ( t_bo_ b(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)457 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 458 458 459 459 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) … … 461 461 dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 462 462 463 q_i_ b(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0)463 q_i_1d(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0) 464 464 465 465 ENDIF ! fc_bo_i … … 476 476 ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K) 477 477 478 zt_i_new = zswitch_sal * t_bo_ b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)478 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 479 479 480 480 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) … … 482 482 & + rcp * ( ztmelts-rtt ) 483 483 484 zEw = rcp * ( t_bo_ b(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)484 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 485 485 486 486 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 487 487 488 488 ! Contribution to heat flux to the ocean [W.m-2], >0 489 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice489 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 490 490 491 491 ! Total heat flux used in this process [W.m-2], <0 492 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_ b(ji) * zdE * r1_rdtice492 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 493 493 494 494 ! Contribution to salt flux, <0 495 sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_ b(ji) * zfmdt * r1_rdtice495 sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_1d(ji) * zfmdt * r1_rdtice 496 496 497 497 ! Contribution to mass flux, <0 498 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_ b(ji) * dh_i_bott(ji) * r1_rdtice498 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice 499 499 500 500 ! update heat content (J.m-2) and layer thickness 501 qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_ b(ji,nlay_i+1)501 qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_1d(ji,nlay_i+1) 502 502 h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji) 503 503 ENDIF … … 512 512 IF( zf_tt(ji) >= 0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared from surface melting 513 513 514 ztmelts = - tmut * s_i_ b(ji,jk) + rtt ! Melting point of layer jk (K)515 516 IF( t_i_ b(ji,jk) >= ztmelts ) THEN !!! Internal melting514 ztmelts = - tmut * s_i_1d(ji,jk) + rtt ! Melting point of layer jk (K) 515 516 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 517 517 zintermelt(ji) = 1._wp 518 518 519 zEi = - q_i_ b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)520 521 !!zEw = rcp * ( t_i_ b(ji,jk) - rtt ) ! Specific enthalpy of meltwater at T = t_i_b(J/kg, <0)519 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 520 521 !!zEw = rcp * ( t_i_1d(ji,jk) - rtt ) ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) 522 522 523 523 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) … … 532 532 533 533 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 534 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_ b(ji) * zEi * r1_rdtice535 536 ! Contribution to salt flux (clem: using sm_i_ b and not s_i_b(jk) is ok)537 sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_ b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice534 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 535 536 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 537 sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 538 538 539 539 ! Contribution to mass flux 540 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice540 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 541 541 542 542 ! update heat content (J.m-2) and layer thickness 543 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_ b(ji,jk)543 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 544 544 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 545 545 546 546 ELSE !!! Basal melting 547 547 548 zEi = - q_i_ b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)548 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 549 549 550 550 zEw = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0) … … 567 567 568 568 ! Contribution to heat flux to the ocean [W.m-2], <0 569 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice570 571 ! Contribution to salt flux (clem: using sm_i_ b and not s_i_b(jk) is ok)572 sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_ b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice569 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 570 571 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 572 sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 573 573 574 574 ! Total heat flux used in this process [W.m-2], >0 575 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_ b(ji) * zdE * r1_rdtice575 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 576 576 577 577 ! Contribution to mass flux 578 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice578 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 579 579 580 580 ! update heat content (J.m-2) and layer thickness 581 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_ b(ji,jk)581 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 582 582 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 583 583 ENDIF … … 602 602 ! 603 603 ! ! excessive energy is sent to lateral ablation 604 ! zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_ b(ji) - epsi20 ) )605 ! zq_1cat(ji) = zinda * rhoic * lfus * at_i_ b(ji) / MAX( 1._wp - at_i_b(ji) , epsi20 ) * zdvres ! J.m-2 >=0604 ! zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_1d(ji) - epsi20 ) ) 605 ! zq_1cat(ji) = zinda * rhoic * lfus * at_i_1d(ji) / MAX( 1._wp - at_i_1d(ji) , epsi20 ) * zdvres ! J.m-2 >=0 606 606 ! 607 607 ! ! correct salt and mass fluxes 608 ! sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_ b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation609 ! wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_ b(ji) * zdvres * r1_rdtice608 ! sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation 609 ! wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdvres * r1_rdtice 610 610 ! ENDIF 611 611 ! END DO … … 616 616 !------------------------------------------- 617 617 DO ji = kideb, kiut 618 ht_i_ b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_bott(ji) )618 ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) ) 619 619 END DO 620 620 … … 627 627 DO ji = kideb, kiut 628 628 zq_rema(ji) = zq_su(ji) + zq_bo(ji) 629 ! zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_ b(ji) ) ) ! =1 if snow629 ! zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) ! =1 if snow 630 630 ! zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zq_s(ji) + epsi20 ) ) 631 631 ! zdeltah (ji,1) = - zindh * zindq * zq_rema(ji) / MAX( zq_s(ji), epsi20 ) 632 ! zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_ b(ji) ) ) ! bound melting632 ! zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 633 633 ! zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1) 634 634 ! dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1) 635 ! ht_s_ b (ji) = ht_s_b(ji) + zdeltah(ji,1)635 ! ht_s_1d (ji) = ht_s_1d(ji) + zdeltah(ji,1) 636 636 ! 637 637 ! zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji) ! update available heat (J.m-2) 638 638 ! ! heat used to melt snow 639 ! hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_ b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0)639 ! hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0) 640 640 ! ! Contribution to mass flux 641 ! wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_ b(ji) * zdeltah(ji,1) * r1_rdtice641 ! wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 642 642 ! 643 643 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 644 644 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 645 hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_ b(ji) ) * r1_rdtice645 hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 646 646 647 647 IF( ln_nicep .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) … … 656 656 DO ji = kideb, kiut 657 657 ! 658 dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_ b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic ) )659 660 ht_i_ b(ji) = ht_i_b(ji) + dh_snowice(ji)661 ht_s_ b(ji) = ht_s_b(ji) - dh_snowice(ji)658 dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic ) ) 659 660 ht_i_1d(ji) = ht_i_1d(ji) + dh_snowice(ji) 661 ht_s_1d(ji) = ht_s_1d(ji) - dh_snowice(ji) 662 662 663 663 ! Salinity of snow ice 664 664 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 665 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_ b(ji)665 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 666 666 667 667 ! entrapment during snow ice formation 668 668 ! new salinity difference stored (to be used in limthd_ent.F90) 669 669 IF ( num_sal == 2 ) THEN 670 zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_ b(ji) - epsi10 ) )670 zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) ) 671 671 ! salinity dif due to snow-ice formation 672 dsm_i_si_1d(ji) = ( zs_snic - sm_i_ b(ji) ) * dh_snowice(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch672 dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi10 ) * zswitch 673 673 ! salinity dif due to bottom growth 674 674 IF ( zf_tt(ji) < 0._wp ) THEN 675 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_ b(ji) ) * dh_i_bott(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch675 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi10 ) * zswitch 676 676 ENDIF 677 677 ENDIF … … 685 685 686 686 ! Contribution to heat flux 687 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice687 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 688 688 689 689 ! Contribution to salt flux 690 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_ b(ji) * zfmdt * r1_rdtice690 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice 691 691 692 692 ! Contribution to mass flux 693 693 ! All snow is thrown in the ocean, and seawater is taken to replace the volume 694 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_ b(ji) * dh_snowice(ji) * rhoic * r1_rdtice695 wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_ b(ji) * dh_snowice(ji) * rhosn * r1_rdtice694 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice 695 wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice 696 696 697 697 ! update heat content (J.m-2) and layer thickness 698 qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_ b(ji,1) + zfmdt * zEw698 qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw 699 699 h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 700 700 701 701 ! Total ablation (to debug) 702 IF( ht_i_ b(ji) <= 0._wp ) a_i_b(ji) = 0._wp702 IF( ht_i_1d(ji) <= 0._wp ) a_i_1d(ji) = 0._wp 703 703 704 704 END DO !ji … … 710 710 !clem bug: we should take snow into account here 711 711 DO ji = kideb, kiut 712 zindh = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_ b(ji) ) )713 t_su_ b(ji) = zindh * t_su_b(ji) + ( 1.0 - zindh ) * rtt712 zindh = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 713 t_su_1d(ji) = zindh * t_su_1d(ji) + ( 1.0 - zindh ) * rtt 714 714 END DO ! ji 715 715 … … 717 717 DO ji = kideb,kiut 718 718 ! mask enthalpy 719 zinda = MAX( 0._wp , SIGN( 1._wp, - ht_s_ b(ji) ) )720 q_s_ b(ji,jk) = ( 1.0 - zinda ) * q_s_b(ji,jk)721 ! recalculate t_s_ b from q_s_b722 t_s_ b(ji,jk) = rtt + ( 1._wp - zinda ) * ( - q_s_b(ji,jk) / ( rhosn * cpic ) + lfus / cpic )719 zinda = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) ) ) 720 q_s_1d(ji,jk) = ( 1.0 - zinda ) * q_s_1d(ji,jk) 721 ! recalculate t_s_1d from q_s_1d 722 t_s_1d(ji,jk) = rtt + ( 1._wp - zinda ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 723 723 END DO 724 724 END DO … … 727 727 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 728 728 CALL wrk_dealloc( jpij, zintermelt ) 729 CALL wrk_dealloc( jpij, jkmax, zdeltah, zh_i )729 CALL wrk_dealloc( jpij, nlay_i+1, zdeltah, zh_i ) 730 730 CALL wrk_dealloc( jpij, icount ) 731 731 ! -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r4901 r4902 75 75 !! 76 76 !! ** Inputs / Ouputs : (global commons) 77 !! surface temperature : t_su_ b78 !! ice/snow temperatures : t_i_ b, t_s_b79 !! ice salinities : s_i_ b77 !! surface temperature : t_su_1d 78 !! ice/snow temperatures : t_i_1d, t_s_1d 79 !! ice salinities : s_i_1d 80 80 !! number of layers in the ice/snow: nlay_i, nlay_s 81 81 !! profile of the ice/snow layers : z_i, z_s 82 !! total ice/snow thickness : ht_i_ b, ht_s_b82 !! total ice/snow thickness : ht_i_1d, ht_s_1d 83 83 !! 84 84 !! ** External : … … 98 98 INTEGER :: ii, ij ! temporary dummy loop index 99 99 INTEGER :: numeq ! current reference number of equation 100 INTEGER :: layer! vertical dummy loop index100 INTEGER :: jk ! vertical dummy loop index 101 101 INTEGER :: nconv ! number of iterations in iterative procedure 102 102 INTEGER :: minnumeqmin, maxnumeqmax … … 114 114 REAL(wp) :: zerritmax ! current maximal error on temperature 115 115 REAL(wp), POINTER, DIMENSION(:) :: ztfs ! ice melting point 116 REAL(wp), POINTER, DIMENSION(:) :: ztsu old! old surface temperature (before the iterative procedure )117 REAL(wp), POINTER, DIMENSION(:) :: ztsu oldit! surface temperature at previous iteration116 REAL(wp), POINTER, DIMENSION(:) :: ztsub ! old surface temperature (before the iterative procedure ) 117 REAL(wp), POINTER, DIMENSION(:) :: ztsubit ! surface temperature at previous iteration 118 118 REAL(wp), POINTER, DIMENSION(:) :: zh_i ! ice layer thickness 119 119 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness … … 129 129 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_i ! Radiation absorbed in the ice 130 130 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_i ! Kappa factor in the ice 131 REAL(wp), POINTER, DIMENSION(:,:) :: zti old! Old temperature in the ice131 REAL(wp), POINTER, DIMENSION(:,:) :: ztib ! Old temperature in the ice 132 132 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_i ! Eta factor in the ice 133 133 REAL(wp), POINTER, DIMENSION(:,:) :: ztitemp ! Temporary temperature in the ice to check the convergence … … 137 137 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_s ! Radiation absorbed in the snow 138 138 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_s ! Kappa factor in the snow 139 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_s 140 REAL(wp), POINTER, DIMENSION(:,:) :: ztstemp 141 REAL(wp), POINTER, DIMENSION(:,:) :: zts old! Temporary temperature in the snow142 REAL(wp), POINTER, DIMENSION(:,:) :: z_s 143 REAL(wp), POINTER, DIMENSION(:,:) :: zindterm ! Independent term144 REAL(wp), POINTER, DIMENSION(:,:) :: zindtbis ! temporary independent term139 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_s ! Eta factor in the snow 140 REAL(wp), POINTER, DIMENSION(:,:) :: ztstemp ! Temporary temperature in the snow to check the convergence 141 REAL(wp), POINTER, DIMENSION(:,:) :: ztsb ! Temporary temperature in the snow 142 REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! Vertical cotes of the layers in the snow 143 REAL(wp), POINTER, DIMENSION(:,:) :: zindterm ! Independent term 144 REAL(wp), POINTER, DIMENSION(:,:) :: zindtbis ! temporary independent term 145 145 REAL(wp), POINTER, DIMENSION(:,:) :: zdiagbis 146 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! tridiagonal system terms146 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! tridiagonal system terms 147 147 ! diag errors on heat 148 148 REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err … … 150 150 ! 151 151 CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow ) 152 CALL wrk_alloc( jpij, ztfs, ztsu old, ztsuoldit, zh_i, zh_s, zfsw )152 CALL wrk_alloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 153 153 CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 154 CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, zti old, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0)155 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, zts old, zeta_s, ztstemp, z_s, kjstart=0)156 CALL wrk_alloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis )157 CALL wrk_alloc( jpij, jkmax+2, 3, ztrid )154 CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0) 155 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0) 156 CALL wrk_alloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis ) 157 CALL wrk_alloc( jpij, nlay_i+3, 3, ztrid ) 158 158 159 159 CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err ) … … 162 162 zdq(:) = 0._wp ; zq_ini(:) = 0._wp 163 163 DO ji = kideb, kiut 164 zq_ini(ji) = ( SUM( q_i_ b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + &165 & SUM( q_s_ b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )164 zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) + & 165 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 166 166 END DO 167 167 … … 172 172 DO ji = kideb , kiut 173 173 ! is there snow or not 174 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_ b(ji) ) ) )174 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) ) ) 175 175 ! surface temperature of fusion 176 176 ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 177 177 ! layer thickness 178 zh_i(ji) = ht_i_ b(ji) / REAL( nlay_i )179 zh_s(ji) = ht_s_ b(ji) / REAL( nlay_s )178 zh_i(ji) = ht_i_1d(ji) / REAL( nlay_i ) 179 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 180 180 END DO 181 181 … … 187 187 z_i(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st ice layer 188 188 189 DO layer= 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer190 DO ji = kideb , kiut 191 z_s(ji, layer) = z_s(ji,layer-1) + ht_s_b(ji) / REAL( nlay_s )192 END DO 193 END DO 194 195 DO layer= 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer196 DO ji = kideb , kiut 197 z_i(ji, layer) = z_i(ji,layer-1) + ht_i_b(ji) / REAL( nlay_i )189 DO jk = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 190 DO ji = kideb , kiut 191 z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) / REAL( nlay_s ) 192 END DO 193 END DO 194 195 DO jk = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 196 DO ji = kideb , kiut 197 z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) / REAL( nlay_i ) 198 198 END DO 199 199 END DO … … 216 216 DO ji = kideb , kiut 217 217 ! switches 218 isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_ b(ji) ) ) )218 isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) ) 219 219 ! hs > 0, isnow = 1 220 220 zhsu (ji) = hnzst ! threshold for the computation of i0 221 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_ b(ji) / zhsu(ji) ) )221 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu(ji) ) ) 222 222 223 223 i0(ji) = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) … … 226 226 ! a function of the cloud cover 227 227 ! 228 !i0(ji) = (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_ b(ji)+10.0)228 !i0(ji) = (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_1d(ji)+10.0) 229 229 !formula used in Cice 230 230 END DO … … 248 248 END DO 249 249 250 DO layer= 1, nlay_s ! Radiation through snow250 DO jk = 1, nlay_s ! Radiation through snow 251 251 DO ji = kideb, kiut 252 252 ! ! radiation transmitted below the layer-th snow layer 253 zradtr_s(ji, layer) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,layer) ) ) )253 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) ) 254 254 ! ! radiation absorbed by the layer-th snow layer 255 zradab_s(ji, layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer)255 zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 256 256 END DO 257 257 END DO … … 261 261 END DO 262 262 263 DO layer= 1, nlay_i ! Radiation through ice263 DO jk = 1, nlay_i ! Radiation through ice 264 264 DO ji = kideb, kiut 265 265 ! ! radiation transmitted below the layer-th ice layer 266 zradtr_i(ji, layer) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,layer) ) ) )266 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) ) 267 267 ! ! radiation absorbed by the layer-th ice layer 268 zradab_i(ji, layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer)268 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 269 269 END DO 270 270 END DO … … 280 280 ! 281 281 DO ji = kideb, kiut ! Old surface temperature 282 ztsu old (ji) = t_su_b(ji) ! temperature at the beg of iter pr.283 ztsu oldit(ji) = t_su_b(ji) ! temperature at the previous iter284 t_su_ b (ji) = MIN( t_su_b(ji), ztfs(ji) - ztsu_err ) ! necessary282 ztsub (ji) = t_su_1d(ji) ! temperature at the beg of iter pr. 283 ztsubit(ji) = t_su_1d(ji) ! temperature at the previous iter 284 t_su_1d (ji) = MIN( t_su_1d(ji), ztfs(ji) - ztsu_err ) ! necessary 285 285 zerrit (ji) = 1000._wp ! initial value of error 286 286 END DO 287 287 288 DO layer= 1, nlay_s ! Old snow temperature289 DO ji = kideb , kiut 290 zts old(ji,layer) = t_s_b(ji,layer)291 END DO 292 END DO 293 294 DO layer= 1, nlay_i ! Old ice temperature295 DO ji = kideb , kiut 296 zti old(ji,layer) = t_i_b(ji,layer)288 DO jk = 1, nlay_s ! Old snow temperature 289 DO ji = kideb , kiut 290 ztsb(ji,jk) = t_s_1d(ji,jk) 291 END DO 292 END DO 293 294 DO jk = 1, nlay_i ! Old ice temperature 295 DO ji = kideb , kiut 296 ztib(ji,jk) = t_i_1d(ji,jk) 297 297 END DO 298 298 END DO … … 311 311 IF( thcon_i_swi == 0 ) THEN ! Untersteiner (1964) formula 312 312 DO ji = kideb , kiut 313 ztcond_i(ji,0) = rcdic + zbeta*s_i_ b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt)313 ztcond_i(ji,0) = rcdic + zbeta*s_i_1d(ji,1) / MIN(-epsi10,t_i_1d(ji,1)-rtt) 314 314 ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin) 315 315 END DO 316 DO layer= 1, nlay_i-1316 DO jk = 1, nlay_i-1 317 317 DO ji = kideb , kiut 318 ztcond_i(ji, layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) / &319 MIN(-2.0_wp * epsi10, t_i_ b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)320 ztcond_i(ji, layer) = MAX(ztcond_i(ji,layer),zkimin)318 ztcond_i(ji,jk) = rcdic + zbeta*( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / & 319 MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt) 320 ztcond_i(ji,jk) = MAX(ztcond_i(ji,jk),zkimin) 321 321 END DO 322 322 END DO … … 325 325 IF( thcon_i_swi == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 326 326 DO ji = kideb , kiut 327 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_ b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt ) &328 & - 0.011_wp * ( t_i_ b(ji,1) - rtt )327 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1)-rtt ) & 328 & - 0.011_wp * ( t_i_1d(ji,1) - rtt ) 329 329 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 330 330 END DO 331 DO layer= 1, nlay_i-1331 DO jk = 1, nlay_i-1 332 332 DO ji = kideb , kiut 333 ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) & 334 & / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) & 335 & - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt ) 336 ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin ) 333 ztcond_i(ji,jk) = rcdic + & 334 & 0.090_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) & 335 & / MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt) & 336 & - 0.0055_wp* ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0*rtt ) 337 ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 337 338 END DO 338 339 END DO 339 340 DO ji = kideb , kiut 340 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_ b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt) &341 & - 0.011_wp * ( t_bo_ b(ji) - rtt )341 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN(-epsi10,t_bo_1d(ji)-rtt) & 342 & - 0.011_wp * ( t_bo_1d(ji) - rtt ) 342 343 ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 343 344 END DO … … 355 356 END DO 356 357 357 DO layer= 1, nlay_s-1358 DO ji = kideb , kiut 359 zkappa_s(ji, layer) = 2.0 * rcdsn / &358 DO jk = 1, nlay_s-1 359 DO ji = kideb , kiut 360 zkappa_s(ji,jk) = 2.0 * rcdsn / & 360 361 MAX(epsi10,2.0*zh_s(ji)) 361 362 END DO 362 363 END DO 363 364 364 DO layer= 1, nlay_i-1365 DO jk = 1, nlay_i-1 365 366 DO ji = kideb , kiut 366 367 !-- Ice kappa factors 367 zkappa_i(ji, layer) = 2.0*ztcond_i(ji,layer)/ &368 zkappa_i(ji,jk) = 2.0*ztcond_i(ji,jk)/ & 368 369 MAX(epsi10,2.0*zh_i(ji)) 369 370 END DO … … 384 385 !------------------------------------------------------------------------------| 385 386 ! 386 DO layer= 1, nlay_i387 DO ji = kideb , kiut 388 ztitemp(ji, layer) = t_i_b(ji,layer)389 zspeche_i(ji, layer) = cpic + zgamma*s_i_b(ji,layer)/ &390 MAX((t_i_ b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10)391 zeta_i(ji, layer) = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), &387 DO jk = 1, nlay_i 388 DO ji = kideb , kiut 389 ztitemp(ji,jk) = t_i_1d(ji,jk) 390 zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ & 391 MAX((t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt),epsi10) 392 zeta_i(ji,jk) = rdt_ice / MAX(rhoic*zspeche_i(ji,jk)*zh_i(ji), & 392 393 epsi10) 393 394 END DO 394 395 END DO 395 396 396 DO layer= 1, nlay_s397 DO ji = kideb , kiut 398 ztstemp(ji, layer) = t_s_b(ji,layer)399 zeta_s(ji, layer) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10)397 DO jk = 1, nlay_s 398 DO ji = kideb , kiut 399 ztstemp(ji,jk) = t_s_1d(ji,jk) 400 zeta_s(ji,jk) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10) 400 401 END DO 401 402 END DO … … 408 409 DO ji = kideb , kiut 409 410 ! update of the non solar flux according to the update in T_su 410 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_ b(ji) - ztsuoldit(ji) )411 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) ) 411 412 END DO 412 413 ENDIF … … 432 433 !!ice interior terms (top equation has the same form as the others) 433 434 434 DO numeq=1, jkmax+2435 DO numeq=1,nlay_i+3 435 436 DO ji = kideb , kiut 436 437 ztrid(ji,numeq,1) = 0. … … 445 446 DO numeq = nlay_s + 2, nlay_s + nlay_i 446 447 DO ji = kideb , kiut 447 layer= numeq - nlay_s - 1448 ztrid(ji,numeq,1) = - zeta_i(ji, layer)*zkappa_i(ji,layer-1)449 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji, layer)*(zkappa_i(ji,layer-1) + &450 zkappa_i(ji, layer))451 ztrid(ji,numeq,3) = - zeta_i(ji, layer)*zkappa_i(ji,layer)452 zindterm(ji,numeq) = zti old(ji,layer) + zeta_i(ji,layer)* &453 zradab_i(ji, layer)448 jk = numeq - nlay_s - 1 449 ztrid(ji,numeq,1) = - zeta_i(ji,jk)*zkappa_i(ji,jk-1) 450 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,jk)*(zkappa_i(ji,jk-1) + & 451 zkappa_i(ji,jk)) 452 ztrid(ji,numeq,3) = - zeta_i(ji,jk)*zkappa_i(ji,jk) 453 zindterm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk)* & 454 zradab_i(ji,jk) 454 455 END DO 455 456 ENDDO … … 462 463 + zkappa_i(ji,nlay_i-1) ) 463 464 ztrid(ji,numeq,3) = 0.0 464 zindterm(ji,numeq) = zti old(ji,nlay_i) + zeta_i(ji,nlay_i)* &465 zindterm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 465 466 ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 466 * t_bo_ b(ji) )467 * t_bo_1d(ji) ) 467 468 ENDDO 468 469 469 470 470 471 DO ji = kideb , kiut 471 IF ( ht_s_ b(ji).gt.0.0 ) THEN472 IF ( ht_s_1d(ji).gt.0.0 ) THEN 472 473 ! 473 474 !------------------------------------------------------------------------------| … … 477 478 !!snow interior terms (bottom equation has the same form as the others) 478 479 DO numeq = 3, nlay_s + 1 479 layer= numeq - 1480 ztrid(ji,numeq,1) = - zeta_s(ji, layer)*zkappa_s(ji,layer-1)481 ztrid(ji,numeq,2) = 1.0 + zeta_s(ji, layer)*( zkappa_s(ji,layer-1) + &482 zkappa_s(ji, layer) )483 ztrid(ji,numeq,3) = - zeta_s(ji, layer)*zkappa_s(ji,layer)484 zindterm(ji,numeq) = zts old(ji,layer) + zeta_s(ji,layer)* &485 zradab_s(ji, layer)480 jk = numeq - 1 481 ztrid(ji,numeq,1) = - zeta_s(ji,jk)*zkappa_s(ji,jk-1) 482 ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,jk)*( zkappa_s(ji,jk-1) + & 483 zkappa_s(ji,jk) ) 484 ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 485 zindterm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk)* & 486 zradab_s(ji,jk) 486 487 END DO 487 488 … … 490 491 ztrid(ji,nlay_s+2,3) = 0.0 491 492 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 492 t_bo_ b(ji)493 t_bo_1d(ji) 493 494 ENDIF 494 495 495 IF ( t_su_ b(ji) .LT. rtt ) THEN496 IF ( t_su_1d(ji) .LT. rtt ) THEN 496 497 497 498 !------------------------------------------------------------------------------| … … 506 507 ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 507 508 ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 508 zindterm(ji,1) = dzf(ji)*t_su_ b(ji) - zf(ji)509 zindterm(ji,1) = dzf(ji)*t_su_1d(ji) - zf(ji) 509 510 510 511 !!first layer of snow equation … … 512 513 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s) 513 514 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 514 zindterm(ji,2) = zts old(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)515 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 515 516 516 517 ELSE … … 529 530 zkappa_s(ji,0) * zg1s ) 530 531 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 531 zindterm(ji,2) = zts old(ji,1) + zeta_s(ji,1) * &532 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * & 532 533 ( zradab_s(ji,1) + & 533 zkappa_s(ji,0) * zg1s * t_su_ b(ji) )534 zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 534 535 ENDIF 535 536 ELSE … … 539 540 !------------------------------------------------------------------------------| 540 541 ! 541 IF (t_su_ b(ji) .LT. rtt) THEN542 IF (t_su_1d(ji) .LT. rtt) THEN 542 543 ! 543 544 !------------------------------------------------------------------------------| … … 553 554 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1 554 555 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 555 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_ b(ji) - zf(ji)556 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji) 556 557 557 558 !!first layer of ice equation … … 560 561 + zkappa_i(ji,0) * zg1 ) 561 562 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1) 562 zindterm(ji,numeqmin(ji)+1)= zti old(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)563 zindterm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1) 563 564 564 565 !!case of only one layer in the ice (surface & ice equations are altered) … … 573 574 ztrid(ji,numeqmin(ji)+1,3) = 0.0 574 575 575 zindterm(ji,numeqmin(ji)+1) = zti old(ji,1) + zeta_i(ji,1)* &576 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_ b(ji) )576 zindterm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1)* & 577 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) ) 577 578 ENDIF 578 579 … … 593 594 zg1) 594 595 ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 595 zindterm(ji,numeqmin(ji)) = zti old(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &596 zkappa_i(ji,0) * zg1 * t_su_ b(ji) )596 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 597 zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 597 598 598 599 !!case of only one layer in the ice (surface & ice equations are altered) … … 602 603 zkappa_i(ji,1)) 603 604 ztrid(ji,numeqmin(ji),3) = 0.0 604 zindterm(ji,numeqmin(ji)) = zti old(ji,1) + zeta_i(ji,1)* &605 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_ b(ji)) &606 + t_su_ b(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0605 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)* & 606 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) & 607 + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0 607 608 ENDIF 608 609 … … 623 624 624 625 maxnumeqmax = 0 625 minnumeqmin = jkmax+4626 minnumeqmin = nlay_i+5 626 627 627 628 DO ji = kideb , kiut … … 632 633 END DO 633 634 634 DO layer= minnumeqmin+1, maxnumeqmax635 DO ji = kideb , kiut 636 numeq = min(max(numeqmin(ji)+1, layer),numeqmax(ji))635 DO jk = minnumeqmin+1, maxnumeqmax 636 DO ji = kideb , kiut 637 numeq = min(max(numeqmin(ji)+1,jk),numeqmax(ji)) 637 638 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 638 639 ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) … … 644 645 DO ji = kideb , kiut 645 646 ! ice temperatures 646 t_i_ b(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))647 t_i_1d(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 647 648 END DO 648 649 649 650 DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1 650 651 DO ji = kideb , kiut 651 layer= numeq - nlay_s - 1652 t_i_ b(ji,layer) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &653 t_i_ b(ji,layer+1))/zdiagbis(ji,numeq)652 jk = numeq - nlay_s - 1 653 t_i_1d(ji,jk) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* & 654 t_i_1d(ji,jk+1))/zdiagbis(ji,numeq) 654 655 END DO 655 656 END DO … … 657 658 DO ji = kideb , kiut 658 659 ! snow temperatures 659 IF (ht_s_ b(ji).GT.0._wp) &660 t_s_ b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &661 * t_i_ b(ji,1))/zdiagbis(ji,nlay_s+1) &662 * MAX(0.0,SIGN(1.0,ht_s_ b(ji)))660 IF (ht_s_1d(ji).GT.0._wp) & 661 t_s_1d(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 662 * t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) & 663 * MAX(0.0,SIGN(1.0,ht_s_1d(ji))) 663 664 664 665 ! surface temperature 665 isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_ b(ji) ) ) )666 ztsu oldit(ji) = t_su_b(ji)667 IF( t_su_ b(ji) < ztfs(ji) ) &668 t_su_ b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1) &669 & + REAL( 1 - isnow(ji) )*t_i_ b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))666 isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_1d(ji) ) ) ) 667 ztsubit(ji) = t_su_1d(ji) 668 IF( t_su_1d(ji) < ztfs(ji) ) & 669 t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1) & 670 & + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 670 671 END DO 671 672 ! … … 677 678 ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd 678 679 DO ji = kideb , kiut 679 t_su_ b(ji) = MAX( MIN( t_su_b(ji) , ztfs(ji) ) , 190._wp )680 zerrit(ji) = ABS( t_su_ b(ji) - ztsuoldit(ji) )681 END DO 682 683 DO layer= 1, nlay_s684 DO ji = kideb , kiut 685 t_s_ b(ji,layer) = MAX( MIN( t_s_b(ji,layer), rtt ), 190._wp )686 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_ b(ji,layer) - ztstemp(ji,layer)))687 END DO 688 END DO 689 690 DO layer= 1, nlay_i691 DO ji = kideb , kiut 692 ztmelt_i = -tmut * s_i_ b(ji,layer) + rtt693 t_i_ b(ji,layer) = MAX(MIN(t_i_b(ji,layer),ztmelt_i), 190._wp)694 zerrit(ji) = MAX(zerrit(ji),ABS(t_i_ b(ji,layer) - ztitemp(ji,layer)))680 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , ztfs(ji) ) , 190._wp ) 681 zerrit(ji) = ABS( t_su_1d(ji) - ztsubit(ji) ) 682 END DO 683 684 DO jk = 1, nlay_s 685 DO ji = kideb , kiut 686 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rtt ), 190._wp ) 687 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_1d(ji,jk) - ztstemp(ji,jk))) 688 END DO 689 END DO 690 691 DO jk = 1, nlay_i 692 DO ji = kideb , kiut 693 ztmelt_i = -tmut * s_i_1d(ji,jk) + rtt 694 t_i_1d(ji,jk) = MAX(MIN(t_i_1d(ji,jk),ztmelt_i), 190._wp) 695 zerrit(ji) = MAX(zerrit(ji),ABS(t_i_1d(ji,jk) - ztitemp(ji,jk))) 695 696 END DO 696 697 END DO … … 717 718 DO ji = kideb, kiut 718 719 ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux) 719 IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_ b(ji) - ztsuold(ji) ) )720 IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) ) 720 721 ! ! surface ice conduction flux 721 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_ b(ji) ) ) )722 fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_ b(ji,1) - t_su_b(ji)) &723 & - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_ b(ji,1) - t_su_b(ji))722 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) ) 723 fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & 724 & - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji)) 724 725 ! ! bottom ice conduction flux 725 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_ b(ji) - t_i_b(ji,nlay_i)) )726 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 726 727 END DO 727 728 … … 730 731 !----------------------------------------- 731 732 DO ji = kideb, kiut 732 IF( t_su_b(ji) < rtt ) THEN ! case T_su < 0degC 733 hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 733 IF( t_su_1d(ji) < rtt ) THEN ! case T_su < 0degC 734 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 735 & ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 734 736 ELSE ! case T_su = 0degC 735 hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 737 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 738 & ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 736 739 ENDIF 737 740 END DO … … 742 745 ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 743 746 DO ji = kideb, kiut 744 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_ b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + &745 & SUM( q_s_ b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )747 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) + & 748 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 746 749 zhfx_err(ji) = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ) 747 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_ b(ji)750 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 748 751 END DO 749 752 … … 767 770 DO ji = kideb, kiut 768 771 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 769 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_ b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) )772 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_1d(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 770 773 END DO 771 774 772 775 ! 773 776 CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 774 CALL wrk_dealloc( jpij, ztfs, ztsu old, ztsuoldit, zh_i, zh_s, zfsw )777 CALL wrk_dealloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 775 778 CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 776 CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 777 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart = 0 ) 778 CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis ) 779 CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid ) 779 CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, & 780 & ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 781 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 782 CALL wrk_dealloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis ) 783 CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid ) 780 784 CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) 781 785 … … 798 802 DO jk = 1, nlay_i ! Sea ice energy of melting 799 803 DO ji = kideb, kiut 800 ztmelts = - tmut * s_i_ b(ji,jk) + rtt801 zindb = MAX( 0._wp , SIGN( 1._wp , -(t_i_ b(ji,jk) - rtt) - epsi10 ) )802 q_i_ b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) ) &803 & + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_ b(ji,jk)-rtt, -epsi10 ) ) &804 ztmelts = - tmut * s_i_1d(ji,jk) + rtt 805 zindb = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi10 ) ) 806 q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) & 807 & + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_1d(ji,jk)-rtt, -epsi10 ) ) & 804 808 & - rcp * ( ztmelts-rtt ) ) 805 809 END DO … … 807 811 DO jk = 1, nlay_s ! Snow energy of melting 808 812 DO ji = kideb, kiut 809 q_s_ b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )813 q_s_1d(ji,jk) = rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) 810 814 END DO 811 815 END DO -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
r4900 r4902 146 146 ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0 147 147 DO ji = kideb, kiut 148 hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_ b(ji) * r1_rdtice * &148 hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_rdtice * & 149 149 & ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( qh_i_old(ji,0:nlay_i+1) ) ) 150 150 END DO -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r4901 r4902 72 72 !! - Computation of variation of ice volume and mass 73 73 !! - Computation of frldb after lateral accretion and 74 !! update ht_s_ b, ht_i_band tbif_1d(:,:)74 !! update ht_s_1d, ht_i_1d and tbif_1d(:,:) 75 75 !!------------------------------------------------------------------------ 76 INTEGER :: ji,jj,jk,jl ,jm! dummy loop indices77 INTEGER :: layer, nbpac! local integers78 INTEGER :: ii, ij, iter ! - -76 INTEGER :: ji,jj,jk,jl ! dummy loop indices 77 INTEGER :: nbpac ! local integers 78 INTEGER :: ii, ij, iter ! - - 79 79 REAL(wp) :: ztmelts, zdv, zfrazb, zweight, zindb, zinda, zde ! local scalars 80 80 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - - … … 90 90 REAL(wp) :: zv_newfra 91 91 92 INTEGER , POINTER, DIMENSION(:) :: jcat ! indexes of categories where new ice grows92 INTEGER , POINTER, DIMENSION(:) :: jcat ! indexes of categories where new ice grows 93 93 REAL(wp), POINTER, DIMENSION(:) :: zswinew ! switch for new ice or not 94 94 … … 102 102 REAL(wp), POINTER, DIMENSION(:) :: zda_res ! residual area in case of excessive heat budget 103 103 REAL(wp), POINTER, DIMENSION(:) :: zat_i_1d ! total ice fraction 104 REAL(wp), POINTER, DIMENSION(:) :: zat_i_lev ! total ice fraction for level ice only (type 1) 105 REAL(wp), POINTER, DIMENSION(:) :: zv_frazb ! accretion of frazil ice at the ice bottom 104 REAL(wp), POINTER, DIMENSION(:) :: zv_frazb ! accretion of frazil ice at the ice bottom 106 105 REAL(wp), POINTER, DIMENSION(:) :: zvrel_1d ! relative ice / frazil velocity (1D vector) 107 106 108 REAL(wp), POINTER, DIMENSION(:,:) :: zv_ old! old volume of ice in category jl109 REAL(wp), POINTER, DIMENSION(:,:) :: za_ old! old area of ice in category jl110 REAL(wp), POINTER, DIMENSION(:,:) :: za_i_1d 111 REAL(wp), POINTER, DIMENSION(:,:) :: zv_i_1d 112 REAL(wp), POINTER, DIMENSION(:,:) :: zoa_i_1d 113 REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_1d 114 115 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_1d 107 REAL(wp), POINTER, DIMENSION(:,:) :: zv_b ! old volume of ice in category jl 108 REAL(wp), POINTER, DIMENSION(:,:) :: za_b ! 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 113 114 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_1d !: 1-D version of e_i 116 115 117 116 REAL(wp), POINTER, DIMENSION(:,:) :: zvrel ! relative ice / frazil velocity … … 120 119 CALL wrk_alloc( jpij, jcat ) ! integer 121 120 CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 122 CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, z at_i_lev, zv_frazb, zvrel_1d )123 CALL wrk_alloc( jpij,jpl, zv_ old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )124 CALL wrk_alloc( jpij, jkmax,jpl, ze_i_1d )121 CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 122 CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 123 CALL wrk_alloc( jpij,nlay_i+1,jpl, ze_i_1d ) 125 124 CALL wrk_alloc( jpi,jpj, zvrel ) 126 125 … … 303 302 304 303 CALL tab_2d_1d( nbpac, qlead_1d (1:nbpac) , qlead , jpi, jpj, npac(1:nbpac) ) 305 CALL tab_2d_1d( nbpac, t_bo_ b(1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) )304 CALL tab_2d_1d( nbpac, t_bo_1d (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) ) 306 305 CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac) , sfx_opw, jpi, jpj, npac(1:nbpac) ) 307 306 CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw, jpi, jpj, npac(1:nbpac) ) 308 CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw, jpi, jpj, npac(1:nbpac) ) 309 CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 307 CALL tab_2d_1d( nbpac, hicol_1d (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 310 308 CALL tab_2d_1d( nbpac, zvrel_1d (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 311 309 … … 320 318 ! Keep old ice areas and volume in memory 321 319 !----------------------------------------- 322 zv_ old(1:nbpac,:) = zv_i_1d(1:nbpac,:)323 za_ old(1:nbpac,:) = za_i_1d(1:nbpac,:)320 zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:) 321 za_b(1:nbpac,:) = za_i_1d(1:nbpac,:) 324 322 !---------------------- 325 323 ! Thickness of new ice … … 328 326 zh_newice(ji) = hiccrit 329 327 END DO 330 IF( fraz_swi == 1 ) zh_newice(1:nbpac) = hicol_ b(1:nbpac)328 IF( fraz_swi == 1 ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 331 329 332 330 !---------------------- … … 352 350 DO ji = 1, nbpac 353 351 ztmelts = - tmut * zs_newice(ji) + rtt ! Melting point (K) 354 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_ b(ji) ) &355 & + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_ b(ji) - rtt, -epsi10 ) ) &352 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_1d(ji) ) & 353 & + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_1d(ji) - rtt, -epsi10 ) ) & 356 354 & - rcp * ( ztmelts - rtt ) ) 357 355 END DO ! ji … … 371 369 zEi = - ze_newice(ji) / rhoic ! specific enthalpy of forming ice [J/kg] 372 370 373 zEw = rcp * ( t_bo_ b(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_b[J/kg]371 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_1d [J/kg] 374 372 ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 375 373 … … 442 440 DO ji = 1, nbpac 443 441 jl = jcat(ji) ! categroy in which new ice is put 444 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_ old(ji,jl) ) ) ! 0 if old ice442 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) ) ! 0 if old ice 445 443 END DO 446 444 … … 450 448 zinda = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 451 449 ze_i_1d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + & 452 & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_ old(ji,jl) ) &450 & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) ) & 453 451 & * zinda / MAX( zv_i_1d(ji,jl), epsi20 ) 454 452 END DO … … 492 490 DO ji = 1, nbpac 493 491 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) ) ! 0 if no ice and 1 if yes 494 zoa_i_1d(ji,jl) = za_ old(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb492 zoa_i_1d(ji,jl) = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb 495 493 END DO 496 494 END DO … … 501 499 DO jl = 1, jpl 502 500 DO ji = 1, nbpac 503 zdv = zv_i_1d(ji,jl) - zv_ old(ji,jl)501 zdv = zv_i_1d(ji,jl) - zv_b(ji,jl) 504 502 zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 505 503 END DO … … 520 518 CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj ) 521 519 CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj ) 522 CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj )523 520 524 521 CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj ) … … 544 541 CALL wrk_dealloc( jpij, jcat ) ! integer 545 542 CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 546 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, z at_i_lev, zv_frazb, zvrel_1d )547 CALL wrk_dealloc( jpij,jpl, zv_ old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )548 CALL wrk_dealloc( jpij, jkmax,jpl, ze_i_1d )543 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 544 CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 545 CALL wrk_dealloc( jpij,nlay_i+1,jpl, ze_i_1d ) 549 546 CALL wrk_dealloc( jpi,jpj, zvrel ) 550 547 ! -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
r4900 r4902 60 60 !--------------------------------------------------------- 61 61 DO ji = kideb, kiut 62 sm_i_ b(ji) = sm_i_b(ji) + dsm_i_se_1d(ji) + dsm_i_si_1d(ji)62 sm_i_1d(ji) = sm_i_1d(ji) + dsm_i_se_1d(ji) + dsm_i_si_1d(ji) 63 63 END DO 64 64 … … 66 66 ! 1) Constant salinity, constant in time | 67 67 !------------------------------------------------------------------------------| 68 !!gm comment: if num_sal = 1 s_i_new, s_i_ b and sm_i_bcan be set to bulk_sal one for all in the initialisation phase !!68 !!gm comment: if num_sal = 1 s_i_new, s_i_1d and sm_i_1d can be set to bulk_sal one for all in the initialisation phase !! 69 69 !!gm ===>>> simplification of almost all test on num_sal value 70 70 IF( num_sal == 1 ) THEN 71 s_i_ b(kideb:kiut,1:nlay_i) = bulk_sal72 sm_i_ b(kideb:kiut) = bulk_sal71 s_i_1d (kideb:kiut,1:nlay_i) = bulk_sal 72 sm_i_1d(kideb:kiut) = bulk_sal 73 73 s_i_new(kideb:kiut) = bulk_sal 74 74 ENDIF … … 83 83 ! Switches 84 84 !---------- 85 iflush = MAX( 0._wp , SIGN( 1._wp , t_su_ b(ji) - rtt ) )! =1 if summer86 igravdr = MAX( 0._wp , SIGN( 1._wp , t_bo_ b(ji) - t_su_b(ji) ) ) ! =1 if t_su < t_bo85 iflush = MAX( 0._wp , SIGN( 1._wp , t_su_1d(ji) - rtt ) ) ! =1 if summer 86 igravdr = MAX( 0._wp , SIGN( 1._wp , t_bo_1d(ji) - t_su_1d(ji) ) ) ! =1 if t_su < t_bo 87 87 88 88 !--------------------- … … 90 90 !--------------------- 91 91 ! drainage by gravity drainage 92 dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_ b(ji) - sal_G , 0._wp ) / time_G * rdt_ice92 dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_1d(ji) - sal_G , 0._wp ) / time_G * rdt_ice 93 93 ! drainage by flushing 94 dsm_i_fl_1d(ji) = - iflush * MAX( sm_i_ b(ji) - sal_F , 0._wp ) / time_F * rdt_ice94 dsm_i_fl_1d(ji) = - iflush * MAX( sm_i_1d(ji) - sal_F , 0._wp ) / time_F * rdt_ice 95 95 96 96 !----------------- … … 99 99 ! only drainage terms ( gravity drainage and flushing ) 100 100 ! snow ice / bottom sources are added in lim_thd_ent to conserve energy 101 sm_i_ b(ji) = sm_i_b(ji) + dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji)101 sm_i_1d(ji) = sm_i_1d(ji) + dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji) 102 102 103 103 !---------------------------- 104 104 ! Salt flux - brine drainage 105 105 !---------------------------- 106 sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoic * a_i_ b(ji) * ht_i_b(ji) * ( dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji) ) * r1_rdtice106 sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoic * a_i_1d(ji) * ht_i_1d(ji) * ( dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji) ) * r1_rdtice 107 107 108 108 END DO -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r4900 r4902 63 63 INTEGER, INTENT(in) :: kt ! number of iteration 64 64 ! 65 INTEGER :: ji, jj, jk, jl, layer! dummy loop indices65 INTEGER :: ji, jj, jk, jl, jn ! dummy loop indices 66 66 INTEGER :: initad ! number of sub-timestep for the advection 67 67 INTEGER :: ierr ! error status … … 85 85 CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold ) 86 86 CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 87 CALL wrk_alloc( jpi, jpj, jkmax, jpl, zs0e )87 CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, zs0e ) 88 88 89 89 CALL wrk_alloc( jpi, jpj, jpl, zaiold, zhimax, zviold, zvsold ) ! clem … … 167 167 168 168 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 169 DO j k= 1,initad169 DO jn = 1,initad 170 170 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area 171 171 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) … … 197 197 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & 198 198 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 199 DO layer= 1, nlay_i !--- ice heat contents ---200 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0e(:,:, layer,jl), sxe (:,:,layer,jl), &201 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &202 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )203 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:, layer,jl), sxe (:,:,layer,jl), &204 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &205 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )199 DO jk = 1, nlay_i !--- ice heat contents --- 200 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), & 201 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 202 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 203 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), & 204 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 205 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 206 206 END DO 207 207 END DO 208 208 END DO 209 209 ELSE 210 DO j k= 1, initad210 DO jn = 1, initad 211 211 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area 212 212 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) … … 239 239 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & 240 240 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 241 DO layer= 1, nlay_i !--- ice heat contents ---242 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0e(:,:, layer,jl), sxe (:,:,layer,jl), &243 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &244 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )245 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:, layer,jl), sxe (:,:,layer,jl), &246 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &247 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )241 DO jk = 1, nlay_i !--- ice heat contents --- 242 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), & 243 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 244 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 245 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), & 246 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 247 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 248 248 END DO 249 249 END DO … … 506 506 CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold ) 507 507 CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 508 CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e )508 CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, zs0e ) 509 509 510 510 CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zaiold, zhimax ) ! clem -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
r4900 r4902 69 69 !! 70 70 !!--------------------------------------------------------------------- 71 INTEGER :: ji, jj, jk, jl, jm ! dummy loop indices 72 INTEGER :: jbnd1, jbnd2 71 INTEGER :: ji, jj, jk, jl ! dummy loop indices 73 72 INTEGER :: i_ice_switch 74 73 REAL(wp) :: zsal … … 93 92 ! Rebin categories with thickness out of bounds 94 93 !---------------------------------------------------- 95 DO jm = 1, jpm 96 jbnd1 = ice_cat_bounds(jm,1) 97 jbnd2 = ice_cat_bounds(jm,2) 98 IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm) 99 END DO 94 IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl) 100 95 101 96 at_i(:,:) = 0._wp … … 126 121 ! Final thickness distribution rebinning 127 122 ! -------------------------------------- 128 DO jm = 1, jpm 129 jbnd1 = ice_cat_bounds(jm,1) 130 jbnd2 = ice_cat_bounds(jm,2) 131 IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm) 132 IF (ice_ncat_types(jm) .EQ. 1 ) THEN 133 ENDIF 134 END DO 123 IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl) 135 124 136 125 !----------------- … … 161 150 ! Diagnostics 162 151 ! ------------------------------------------------- 163 d_u_ice_dyn(:,:) = u_ice(:,:) - old_u_ice(:,:)164 d_v_ice_dyn(:,:) = v_ice(:,:) - old_v_ice(:,:)165 d_a_i_trp (:,:,:) = a_i (:,:,:) - old_a_i(:,:,:)166 d_v_s_trp (:,:,:) = v_s (:,:,:) - old_v_s(:,:,:)167 d_v_i_trp (:,:,:) = v_i (:,:,:) - old_v_i(:,:,:)168 d_e_s_trp (:,:,:,:) = e_s (:,:,:,:) - old_e_s(:,:,:,:)169 d_e_i_trp (:,:,1:nlay_i,:) = e_i (:,:,1:nlay_i,:) - old_e_i(:,:,1:nlay_i,:)170 d_oa_i_trp (:,:,:) = oa_i (:,:,:) - o ld_oa_i(:,:,:)152 d_u_ice_dyn(:,:) = u_ice(:,:) - u_ice_b(:,:) 153 d_v_ice_dyn(:,:) = v_ice(:,:) - v_ice_b(:,:) 154 d_a_i_trp (:,:,:) = a_i (:,:,:) - a_i_b (:,:,:) 155 d_v_s_trp (:,:,:) = v_s (:,:,:) - v_s_b (:,:,:) 156 d_v_i_trp (:,:,:) = v_i (:,:,:) - v_i_b (:,:,:) 157 d_e_s_trp (:,:,:,:) = e_s (:,:,:,:) - e_s_b (:,:,:,:) 158 d_e_i_trp (:,:,1:nlay_i,:) = e_i (:,:,1:nlay_i,:) - e_i_b(:,:,1:nlay_i,:) 159 d_oa_i_trp (:,:,:) = oa_i (:,:,:) - oa_i_b (:,:,:) 171 160 d_smv_i_trp(:,:,:) = 0._wp 172 IF( num_sal == 2 ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)161 IF( num_sal == 2 ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - smv_i_b(:,:,:) 173 162 174 163 ! conservation test … … 186 175 CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_update1 : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :') 187 176 CALL prt_ctl(tab2d_1=d_u_ice_dyn, clinfo1=' lim_update1 : d_u_ice_dyn :', tab2d_2=d_v_ice_dyn, clinfo2=' d_v_ice_dyn :') 188 CALL prt_ctl(tab2d_1= old_u_ice , clinfo1=' lim_update1 : old_u_ice :', tab2d_2=old_v_ice , clinfo2=' old_v_ice:')177 CALL prt_ctl(tab2d_1=u_ice_b , clinfo1=' lim_update1 : u_ice_b :', tab2d_2=v_ice_b , clinfo2=' v_ice_b :') 189 178 190 179 DO jl = 1, jpl … … 199 188 CALL prt_ctl(tab2d_1=o_i (:,:,jl) , clinfo1= ' lim_update1 : o_i : ') 200 189 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_update1 : a_i : ') 201 CALL prt_ctl(tab2d_1= old_a_i (:,:,jl) , clinfo1= ' lim_update1 : old_a_i: ')190 CALL prt_ctl(tab2d_1=a_i_b (:,:,jl) , clinfo1= ' lim_update1 : a_i_b : ') 202 191 CALL prt_ctl(tab2d_1=d_a_i_trp (:,:,jl) , clinfo1= ' lim_update1 : d_a_i_trp : ') 203 192 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_update1 : v_i : ') 204 CALL prt_ctl(tab2d_1= old_v_i (:,:,jl) , clinfo1= ' lim_update1 : old_v_i: ')193 CALL prt_ctl(tab2d_1=v_i_b (:,:,jl) , clinfo1= ' lim_update1 : v_i_b : ') 205 194 CALL prt_ctl(tab2d_1=d_v_i_trp (:,:,jl) , clinfo1= ' lim_update1 : d_v_i_trp : ') 206 195 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_update1 : v_s : ') 207 CALL prt_ctl(tab2d_1= old_v_s (:,:,jl) , clinfo1= ' lim_update1 : old_v_s: ')196 CALL prt_ctl(tab2d_1=v_s_b (:,:,jl) , clinfo1= ' lim_update1 : v_s_b : ') 208 197 CALL prt_ctl(tab2d_1=d_v_s_trp (:,:,jl) , clinfo1= ' lim_update1 : d_v_s_trp : ') 209 198 CALL prt_ctl(tab2d_1=e_i (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1 : e_i1 : ') 210 CALL prt_ctl(tab2d_1= old_e_i (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1 : old_e_i1: ')199 CALL prt_ctl(tab2d_1=e_i_b (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1 : e_i1_b : ') 211 200 CALL prt_ctl(tab2d_1=d_e_i_trp (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1 : de_i1_trp : ') 212 201 CALL prt_ctl(tab2d_1=e_i (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1 : e_i2 : ') 213 CALL prt_ctl(tab2d_1= old_e_i (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1 : old_e_i2: ')202 CALL prt_ctl(tab2d_1=e_i_b (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1 : e_i2_b : ') 214 203 CALL prt_ctl(tab2d_1=d_e_i_trp (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1 : de_i2_trp : ') 215 204 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_update1 : e_snow : ') 216 CALL prt_ctl(tab2d_1= old_e_s (:,:,1,jl) , clinfo1= ' lim_update1 : old_e_snow: ')205 CALL prt_ctl(tab2d_1=e_s_b (:,:,1,jl) , clinfo1= ' lim_update1 : e_snow_b : ') 217 206 CALL prt_ctl(tab2d_1=d_e_s_trp (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1 : d_e_s_trp : ') 218 207 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_update1 : smv_i : ') 219 CALL prt_ctl(tab2d_1= old_smv_i (:,:,jl) , clinfo1= ' lim_update1 : old_smv_i: ')208 CALL prt_ctl(tab2d_1=smv_i_b (:,:,jl) , clinfo1= ' lim_update1 : smv_i_b : ') 220 209 CALL prt_ctl(tab2d_1=d_smv_i_trp(:,:,jl) , clinfo1= ' lim_update1 : d_smv_i_trp : ') 221 210 CALL prt_ctl(tab2d_1=oa_i (:,:,jl) , clinfo1= ' lim_update1 : oa_i : ') 222 CALL prt_ctl(tab2d_1=o ld_oa_i (:,:,jl) , clinfo1= ' lim_update1 : old_oa_i: ')211 CALL prt_ctl(tab2d_1=oa_i_b (:,:,jl) , clinfo1= ' lim_update1 : oa_i_b : ') 223 212 CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl) , clinfo1= ' lim_update1 : d_oa_i_trp : ') 224 213 -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
r4900 r4902 67 67 !! 68 68 !!--------------------------------------------------------------------- 69 INTEGER :: ji, jj, jk, jl, jm ! dummy loop indices 70 INTEGER :: jbnd1, jbnd2 69 INTEGER :: ji, jj, jk, jl ! dummy loop indices 71 70 INTEGER :: i_ice_switch 72 71 REAL(wp) :: zh, zsal … … 89 88 ! Rebin categories with thickness out of bounds 90 89 !---------------------------------------------------- 91 DO jm = 1, jpm 92 jbnd1 = ice_cat_bounds(jm,1) 93 jbnd2 = ice_cat_bounds(jm,2) 94 IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm) 95 END DO 90 IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl) 96 91 97 92 !---------------------------------------------------------------------- 98 93 ! Constrain the thickness of the smallest category above hiclim 99 94 !---------------------------------------------------------------------- 100 DO jm = 1, jpm 101 DO jj = 1, jpj 102 DO ji = 1, jpi 103 jl = ice_cat_bounds(jm,1) 104 IF( v_i(ji,jj,jl) > 0._wp .AND. ht_i(ji,jj,jl) < hiclim ) THEN 105 zh = hiclim / ht_i(ji,jj,jl) 106 ht_s(ji,jj,jl) = ht_s(ji,jj,jl) * zh 107 ht_i(ji,jj,jl) = ht_i(ji,jj,jl) * zh 108 a_i (ji,jj,jl) = a_i(ji,jj,jl) / zh 109 ENDIF 110 END DO !ji 111 END DO !jj 112 END DO !jm 95 DO jj = 1, jpj 96 DO ji = 1, jpi 97 IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < hiclim ) THEN 98 zh = hiclim / ht_i(ji,jj,1) 99 ht_s(ji,jj,1) = ht_s(ji,jj,1) * zh 100 ht_i(ji,jj,1) = ht_i(ji,jj,1) * zh 101 a_i (ji,jj,1) = a_i(ji,jj,1) / zh 102 ENDIF 103 END DO 104 END DO 113 105 114 106 !----------------------------------------------------- … … 139 131 ! Final thickness distribution rebinning 140 132 ! -------------------------------------- 141 DO jm = 1, jpm 142 jbnd1 = ice_cat_bounds(jm,1) 143 jbnd2 = ice_cat_bounds(jm,2) 144 IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm) 145 IF (ice_ncat_types(jm) .EQ. 1 ) THEN 146 ENDIF 147 END DO 133 IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl ) 148 134 149 135 !----------------- … … 196 182 ! Diagnostics 197 183 ! ------------------------------------------------- 198 d_a_i_thd(:,:,:) = a_i(:,:,:) - old_a_i(:,:,:)199 d_v_s_thd(:,:,:) = v_s(:,:,:) - old_v_s(:,:,:)200 d_v_i_thd(:,:,:) = v_i(:,:,:) - old_v_i(:,:,:)201 d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:)202 d_e_i_thd(:,:,1:nlay_i,:) = e_i(:,:,1:nlay_i,:) - old_e_i(:,:,1:nlay_i,:)203 !?? d_oa_i_thd(:,:,:) = oa_i (:,:,:) - o ld_oa_i(:,:,:)184 d_a_i_thd(:,:,:) = a_i(:,:,:) - a_i_b(:,:,:) 185 d_v_s_thd(:,:,:) = v_s(:,:,:) - v_s_b(:,:,:) 186 d_v_i_thd(:,:,:) = v_i(:,:,:) - v_i_b(:,:,:) 187 d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - e_s_b(:,:,:,:) 188 d_e_i_thd(:,:,1:nlay_i,:) = e_i(:,:,1:nlay_i,:) - e_i_b(:,:,1:nlay_i,:) 189 !?? d_oa_i_thd(:,:,:) = oa_i (:,:,:) - oa_i_b (:,:,:) 204 190 d_smv_i_thd(:,:,:) = 0._wp 205 IF( num_sal == 2 ) d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)191 IF( num_sal == 2 ) d_smv_i_thd(:,:,:) = smv_i(:,:,:) - smv_i_b(:,:,:) 206 192 ! diag only (clem) 207 193 dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday … … 211 197 DO ji = 1, jpi 212 198 diag_heat_dhc(ji,jj) = ( SUM( d_e_i_trp(ji,jj,1:nlay_i,:) + d_e_i_thd(ji,jj,1:nlay_i,:) ) + & 213 & SUM( d_e_s_trp(ji,jj,1:nlay_s,:) + d_e_s_thd(ji,jj,1:nlay_s,:) ) ) * unit_fac * r1_rdtice / area(ji,jj) 199 & SUM( d_e_s_trp(ji,jj,1:nlay_s,:) + d_e_s_thd(ji,jj,1:nlay_s,:) ) & 200 & ) * unit_fac * r1_rdtice / area(ji,jj) 214 201 END DO 215 202 END DO … … 228 215 CALL prt_ctl(tab2d_1=strength , clinfo1=' lim_update2 : strength :') 229 216 CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_update2 : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :') 230 CALL prt_ctl(tab2d_1= old_u_ice , clinfo1=' lim_update2 : old_u_ice :', tab2d_2=old_v_ice , clinfo2=' old_v_ice:')217 CALL prt_ctl(tab2d_1=u_ice_b , clinfo1=' lim_update2 : u_ice_b :', tab2d_2=v_ice_b , clinfo2=' v_ice_b :') 231 218 232 219 DO jl = 1, jpl … … 241 228 CALL prt_ctl(tab2d_1=o_i (:,:,jl) , clinfo1= ' lim_update2 : o_i : ') 242 229 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_update2 : a_i : ') 243 CALL prt_ctl(tab2d_1= old_a_i (:,:,jl) , clinfo1= ' lim_update2 : old_a_i: ')230 CALL prt_ctl(tab2d_1=a_i_b (:,:,jl) , clinfo1= ' lim_update2 : a_i_b : ') 244 231 CALL prt_ctl(tab2d_1=d_a_i_thd (:,:,jl) , clinfo1= ' lim_update2 : d_a_i_thd : ') 245 232 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_update2 : v_i : ') 246 CALL prt_ctl(tab2d_1= old_v_i (:,:,jl) , clinfo1= ' lim_update2 : old_v_i: ')233 CALL prt_ctl(tab2d_1=v_i_b (:,:,jl) , clinfo1= ' lim_update2 : v_i_b : ') 247 234 CALL prt_ctl(tab2d_1=d_v_i_thd (:,:,jl) , clinfo1= ' lim_update2 : d_v_i_thd : ') 248 235 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_update2 : v_s : ') 249 CALL prt_ctl(tab2d_1= old_v_s (:,:,jl) , clinfo1= ' lim_update2 : old_v_s: ')236 CALL prt_ctl(tab2d_1=v_s_b (:,:,jl) , clinfo1= ' lim_update2 : v_s_b : ') 250 237 CALL prt_ctl(tab2d_1=d_v_s_thd (:,:,jl) , clinfo1= ' lim_update2 : d_v_s_thd : ') 251 238 CALL prt_ctl(tab2d_1=e_i (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2 : e_i1 : ') 252 CALL prt_ctl(tab2d_1= old_e_i (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2 : old_e_i1: ')239 CALL prt_ctl(tab2d_1=e_i_b (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2 : e_i1_b : ') 253 240 CALL prt_ctl(tab2d_1=d_e_i_thd (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2 : de_i1_thd : ') 254 241 CALL prt_ctl(tab2d_1=e_i (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2 : e_i2 : ') 255 CALL prt_ctl(tab2d_1= old_e_i (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2 : old_e_i2: ')242 CALL prt_ctl(tab2d_1=e_i_b (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2 : e_i2_b : ') 256 243 CALL prt_ctl(tab2d_1=d_e_i_thd (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2 : de_i2_thd : ') 257 244 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_update2 : e_snow : ') 258 CALL prt_ctl(tab2d_1= old_e_s (:,:,1,jl) , clinfo1= ' lim_update2 : old_e_snow: ')245 CALL prt_ctl(tab2d_1=e_s_b (:,:,1,jl) , clinfo1= ' lim_update2 : e_snow_b : ') 259 246 CALL prt_ctl(tab2d_1=d_e_s_thd (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2 : d_e_s_thd : ') 260 247 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_update2 : smv_i : ') 261 CALL prt_ctl(tab2d_1= old_smv_i (:,:,jl) , clinfo1= ' lim_update2 : old_smv_i: ')248 CALL prt_ctl(tab2d_1=smv_i_b (:,:,jl) , clinfo1= ' lim_update2 : smv_i_b : ') 262 249 CALL prt_ctl(tab2d_1=d_smv_i_thd(:,:,jl) , clinfo1= ' lim_update2 : d_smv_i_thd : ') 263 250 CALL prt_ctl(tab2d_1=oa_i (:,:,jl) , clinfo1= ' lim_update2 : oa_i : ') 264 CALL prt_ctl(tab2d_1=o ld_oa_i (:,:,jl) , clinfo1= ' lim_update2 : old_oa_i: ')251 CALL prt_ctl(tab2d_1=oa_i_b (:,:,jl) , clinfo1= ' lim_update2 : oa_i_b : ') 265 252 CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl) , clinfo1= ' lim_update2 : d_oa_i_thd : ') 266 253 -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r4900 r4902 464 464 ! Vertically constant, constant in time 465 465 !--------------------------------------- 466 IF( num_sal == 1 ) s_i_ b(:,:) = bulk_sal466 IF( num_sal == 1 ) s_i_1d(:,:) = bulk_sal 467 467 468 468 !------------------------------------------------------ … … 473 473 ! 474 474 DO ji = kideb, kiut ! Slope of the linear profile zs_zero 475 z_slope_s(ji) = 2._wp * sm_i_ b(ji) / MAX( epsi10 , ht_i_b(ji) )475 z_slope_s(ji) = 2._wp * sm_i_1d(ji) / MAX( epsi10 , ht_i_1d(ji) ) 476 476 END DO 477 477 … … 489 489 ij = ( npb(ji) - 1 ) / jpi + 1 490 490 ! zind0 = 1 if sm_i le s_i_0 and 0 otherwise 491 zind0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i_ b(ji) ) )491 zind0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i_1d(ji) ) ) 492 492 ! zind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws 493 zind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_ b(ji) ) )493 zind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_1d(ji) ) ) 494 494 ! if 2.sm_i GE sss_m then zindbal = 1 495 495 ! this is to force a constant salinity profile in the Baltic Sea 496 zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_ b(ji) - sss_m(ii,ij) ) )496 zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) ) 497 497 ! 498 zalpha = ( zind0 + zind01 * ( sm_i_ b(ji) * dummy_fac0 + dummy_fac1 ) ) * ( 1.0 - zindbal )498 zalpha = ( zind0 + zind01 * ( sm_i_1d(ji) * dummy_fac0 + dummy_fac1 ) ) * ( 1.0 - zindbal ) 499 499 ! 500 zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_ b(ji) * dummy_fac2500 zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * dummy_fac2 501 501 ! weighting the profile 502 s_i_ b(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_b(ji)502 s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji) 503 503 END DO ! ji 504 504 END DO ! jk … … 512 512 IF( num_sal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 513 513 ! 514 sm_i_ b(:) = 2.30_wp514 sm_i_1d(:) = 2.30_wp 515 515 ! 516 516 !CDIR NOVERRCHK … … 519 519 zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) ) 520 520 DO ji = kideb, kiut 521 s_i_ b(ji,jk) = zsal521 s_i_1d(ji,jk) = zsal 522 522 END DO 523 523 END DO -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
r4900 r4902 298 298 !!---------------------------------------------------------------------- 299 299 300 CALL histdef( kid, "iicethic", "Ice thickness" , "m" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 301 CALL histdef( kid, "iiceconc", "Ice concentration" , "%" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 302 CALL histdef( kid, "iicetemp", "Ice temperature" , "C" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 303 CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 304 CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 305 CALL histdef( kid, "iicestru", "i-Wind stress over ice (I-pt)", "Pa", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 306 CALL histdef( kid, "iicestrv", "j-Wind stress over ice (I-pt)", "Pa", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 307 CALL histdef( kid, "iicesflx", "Solar flux over ocean" , "w/m2" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 308 CALL histdef( kid, "iicenflx", "Non-solar flux over ocean" , "w/m2" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 309 CALL histdef( kid, "isnowpre", "Snow precipitation" , "kg/m2/s", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 310 CALL histdef( kid, "iicesali", "Ice salinity" , "PSU" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 311 CALL histdef( kid, "iicevolu", "Ice volume" , "m" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 312 CALL histdef( kid, "iicedive", "Ice divergence" , "10-8s-1", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 313 CALL histdef( kid, "iicebopr", "Ice bottom production" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 314 CALL histdef( kid, "iicedypr", "Ice dynamic production" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 315 CALL histdef( kid, "iicelapr", "Ice open water prod" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 316 CALL histdef( kid, "iicesipr", "Snow ice production " , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 317 CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 318 CALL histdef( kid, "iicebome", "Ice bottom melt" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 319 CALL histdef( kid, "iicesume", "Ice surface melt" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 320 CALL histdef( kid, "iisfxdyn", "Salt flux from dynmics" , "" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 321 CALL histdef( kid, "iisfxres", "Salt flux from limupdate", "" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 300 CALL histdef( kid, "iicethic", "Ice thickness" , "m" , & 301 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 302 CALL histdef( kid, "iiceconc", "Ice concentration" , "%" , & 303 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 304 CALL histdef( kid, "iicetemp", "Ice temperature" , "C" , & 305 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 306 CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)" , "m/s" , & 307 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 308 CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)" , "m/s" , & 309 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 310 CALL histdef( kid, "iicestru", "i-Wind stress over ice (I-pt)", "Pa", & 311 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 312 CALL histdef( kid, "iicestrv", "j-Wind stress over ice (I-pt)", "Pa", & 313 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 314 CALL histdef( kid, "iicesflx", "Solar flux over ocean" , "w/m2" , & 315 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 316 CALL histdef( kid, "iicenflx", "Non-solar flux over ocean" , "w/m2" , & 317 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 318 CALL histdef( kid, "isnowpre", "Snow precipitation" , "kg/m2/s", & 319 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 320 CALL histdef( kid, "iicesali", "Ice salinity" , "PSU" , & 321 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 322 CALL histdef( kid, "iicevolu", "Ice volume" , "m" , & 323 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 324 CALL histdef( kid, "iicedive", "Ice divergence" , "10-8s-1", & 325 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 326 CALL histdef( kid, "iicebopr", "Ice bottom production" , "m/s" , & 327 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 328 CALL histdef( kid, "iicedypr", "Ice dynamic production" , "m/s" , & 329 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 330 CALL histdef( kid, "iicelapr", "Ice open water prod" , "m/s" , & 331 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 332 CALL histdef( kid, "iicesipr", "Snow ice production " , "m/s" , & 333 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 334 CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s" , & 335 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 336 CALL histdef( kid, "iicebome", "Ice bottom melt" , "m/s" , & 337 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 338 CALL histdef( kid, "iicesume", "Ice surface melt" , "m/s" , & 339 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 340 CALL histdef( kid, "iisfxdyn", "Salt flux from dynmics" , "" , & 341 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 342 CALL histdef( kid, "iisfxres", "Salt flux from limupdate", "" , & 343 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 322 344 323 345 CALL histend( kid, snc4set ) ! end of the file definition -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90
r4900 r4902 12 12 13 13 ! !!! ice thermodynamics 14 INTEGER, PUBLIC, PARAMETER :: jkmax = 6 !: maximum number of ice layers15 14 INTEGER, PUBLIC, PARAMETER :: nlay_i = 5 !: number of ice layers 16 15 INTEGER, PUBLIC, PARAMETER :: nlay_s = 1 !: number of snow layers … … 18 17 ! !!! ice mechanical redistribution 19 18 INTEGER, PUBLIC, PARAMETER :: jpl = 5 !: number of ice categories 20 INTEGER, PUBLIC, PARAMETER :: jpm = 1 !: number of ice types21 19 22 20 !!---------------------------------------------------------------------- -
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
r4901 r4902 34 34 !!----------------------------- 35 35 !: In ice thermodynamics, to spare memory, the vectors are folded 36 !: from 1D to 2D vectors. The following variables, with ending _1d (or _b)36 !: from 1D to 2D vectors. The following variables, with ending _1d 37 37 !: are the variables corresponding to 2d vectors 38 38 … … 40 40 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: npac !: correspondance between points (lateral accretion) 41 41 42 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qlead_1d !: <==> the 2D qlead43 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ftr_ice_1d !: <==> the 2D ftr_ice44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qsr_ice_1d !: <==> the 2D qsr_ice45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fr1_i0_1d !: <==> the 2D fr1_i046 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fr2_i0_1d !: <==> the 2D fr2_i047 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qns_ice_1d !: <==> the 2D qns_ice48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: t_bo_ b !: <==> the 2D t_bo42 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qlead_1d 43 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ftr_ice_1d 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qsr_ice_1d 45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fr1_i0_1d 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fr2_i0_1d 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qns_ice_1d 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: t_bo_1d 49 49 50 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_sum_1d … … 65 65 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_res_1d 66 66 67 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_ice_1d !: <==> the 2D wfx_ice 68 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_snw_1d !: <==> the 2D wfx_snw 69 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_sub_1d !: <==> the 2D wfx_sub 67 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_snw_1d 68 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_sub_1d 70 69 71 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_bog_1d !: <==> the 2D wfx_ice72 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_bom_1d !: <==> the 2D wfx_ice73 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_sum_1d !: <==> the 2D wfx_ice74 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_sni_1d !: <==> the 2D wfx_ice75 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_opw_1d !: <==> the 2D wfx_ice76 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_res_1d !: <==> the 2D wfx_ice77 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_spr_1d !: <==> the 2D wfx_ice70 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_bog_1d 71 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_bom_1d 72 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_sum_1d 73 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_sni_1d 74 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_opw_1d 75 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_res_1d 76 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_spr_1d 78 77 79 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_bri_1d !: <==> the 2D sfx_bri80 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_bog_1d !:81 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_bom_1d !:82 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_sum_1d !:83 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_sni_1d !:84 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_opw_1d !:85 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_res_1d !:78 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_bri_1d 79 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_bog_1d 80 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_bom_1d 81 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_sum_1d 82 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_sni_1d 83 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_opw_1d 84 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_res_1d 86 85 87 86 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sprecip_1d !: <==> the 2D sprecip 88 87 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: frld_1d !: <==> the 2D frld 89 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: at_i_ b!: <==> the 2D at_i90 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fhtur_1d 88 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: at_i_1d !: <==> the 2D at_i 89 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fhtur_1d !: <==> the 2D fhtur 91 90 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fhld_1d !: <==> the 2D fhld 92 91 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dqns_ice_1d !: <==> the 2D dqns_ice … … 100 99 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dsm_i_se_1d !: Ice salinity variations due to basal salt entrapment 101 100 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dsm_i_si_1d !: Ice salinity variations due to lateral accretion 102 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hicol_ b !: Ice collection thickness accumulated in fleads101 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hicol_1d !: Ice collection thickness accumulated in leads 103 102 104 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: t_su_ b!: <==> the 2D t_su105 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_i_ b!: <==> the 2D a_i106 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_i_ b!: <==> the 2D ht_s107 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_s_ b!: <==> the 2D ht_i108 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fc_su !: Surface Conduction flux109 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fc_bo_i !: Bottom Conduction flux110 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_s_tot !: Snow accretion/ablation [m]111 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_i_surf !: Ice surface accretion/ablation [m]112 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_i_bott !: Ice bottom accretion/ablation [m]113 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_snowice !: Snow ice formation [m of ice]114 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sm_i_ b!: Ice bulk salinity [ppt]115 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: s_i_new !: Salinity of new ice at the bottom103 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: t_su_1d !: <==> the 2D t_su 104 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_i_1d !: <==> the 2D a_i 105 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_i_1d !: <==> the 2D ht_s 106 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_s_1d !: <==> the 2D ht_i 107 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fc_su !: Surface Conduction flux 108 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fc_bo_i !: Bottom Conduction flux 109 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_s_tot !: Snow accretion/ablation [m] 110 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_i_surf !: Ice surface accretion/ablation [m] 111 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_i_bott !: Ice bottom accretion/ablation [m] 112 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_snowice !: Snow ice formation [m of ice] 113 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sm_i_1d !: Ice bulk salinity [ppt] 114 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: s_i_new !: Salinity of new ice at the bottom 116 115 117 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_s_ b!: corresponding to the 2D var t_s118 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_i_ b!: corresponding to the 2D var t_i119 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: s_i_ b!: profiled ice salinity120 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_i_ b!: Ice enthalpy per unit volume121 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_s_ b!: Snow enthalpy per unit volume116 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_s_1d !: corresponding to the 2D var t_s 117 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_i_1d !: corresponding to the 2D var t_i 118 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: s_i_1d !: profiled ice salinity 119 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_i_1d !: Ice enthalpy per unit volume 120 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_s_1d !: Snow enthalpy per unit volume 122 121 123 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qh_i_old 124 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_i_old 122 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qh_i_old !: ice heat content (q*h, J.m-2) 123 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_i_old !: ice thickness layer (m) 125 124 126 125 INTEGER , PUBLIC :: jiindex_1d ! 1D index of debugging point … … 146 145 & qsr_ice_1d (jpij) , & 147 146 & fr1_i0_1d(jpij) , fr2_i0_1d(jpij) , qns_ice_1d(jpij) , & 148 & t_bo_b (jpij) , & 149 & hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) ,hfx_dif_1d(jpij) ,hfx_opw_1d(jpij) , & 147 & t_bo_1d (jpij) , & 148 & hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) , & 149 & hfx_dif_1d(jpij) ,hfx_opw_1d(jpij) , & 150 150 & hfx_thd_1d(jpij) , hfx_spr_1d(jpij) , & 151 & hfx_snw_1d(jpij) , hfx_sub_1d(jpij) , hfx_err_1d(jpij) , hfx_res_1d(jpij) , hfx_err_rem_1d(jpij), STAT=ierr(1) ) 151 & hfx_snw_1d(jpij) , hfx_sub_1d(jpij) , hfx_err_1d(jpij) , & 152 & hfx_res_1d(jpij) , hfx_err_rem_1d(jpij), STAT=ierr(1) ) 152 153 ! 153 ALLOCATE( sprecip_1d (jpij) , frld_1d (jpij) , at_i_b (jpij) , & 154 & fhtur_1d (jpij) , wfx_ice_1d (jpij) , wfx_snw_1d (jpij) , wfx_spr_1d (jpij) , & 155 & fhld_1d (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d(jpij) , wfx_bom_1d(jpij) , wfx_sum_1d(jpij) , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d (jpij) , & 154 ALLOCATE( sprecip_1d (jpij) , frld_1d (jpij) , at_i_1d (jpij) , & 155 & fhtur_1d (jpij) , wfx_snw_1d (jpij) , wfx_spr_1d (jpij) , & 156 & fhld_1d (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d(jpij) , wfx_bom_1d(jpij) , & 157 & wfx_sum_1d(jpij) , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d (jpij) , & 156 158 & dqns_ice_1d(jpij) , qla_ice_1d (jpij) , dqla_ice_1d(jpij) , & 157 159 & tatm_ice_1d(jpij) , & 158 160 & i0 (jpij) , & 159 & sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) ,sfx_sum_1d (jpij) ,sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , & 161 & sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) ,sfx_sum_1d (jpij) , & 162 & sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , & 160 163 & dsm_i_fl_1d(jpij) , dsm_i_gd_1d(jpij) , dsm_i_se_1d(jpij) , & 161 & dsm_i_si_1d(jpij) , hicol_ b(jpij) , STAT=ierr(2) )164 & dsm_i_si_1d(jpij) , hicol_1d (jpij) , STAT=ierr(2) ) 162 165 ! 163 ALLOCATE( t_su_ b (jpij) , a_i_b (jpij) , ht_i_b(jpij) , &164 & ht_s_ b(jpij) , fc_su (jpij) , fc_bo_i (jpij) , &166 ALLOCATE( t_su_1d (jpij) , a_i_1d (jpij) , ht_i_1d (jpij) , & 167 & ht_s_1d (jpij) , fc_su (jpij) , fc_bo_i (jpij) , & 165 168 & dh_s_tot (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) , & 166 & dh_snowice(jpij) , sm_i_ b(jpij) , s_i_new (jpij) , &167 & t_s_ b(jpij,nlay_s), &168 & t_i_ b(jpij,jkmax), s_i_b(jpij,jkmax) , &169 & q_i_ b(jpij,jkmax), q_s_b(jpij,jkmax) , &170 & qh_i_old(jpij,0: jkmax), h_i_old(jpij,0:jkmax) , STAT=ierr(3))169 & dh_snowice(jpij) , sm_i_1d (jpij) , s_i_new (jpij) , & 170 & t_s_1d(jpij,nlay_s), & 171 & t_i_1d(jpij,nlay_i+1), s_i_1d(jpij,nlay_i+1) , & 172 & q_i_1d(jpij,nlay_i+1), q_s_1d(jpij,nlay_i+1) , & 173 & qh_i_old(jpij,0:nlay_i+1), h_i_old(jpij,0:nlay_i+1) , STAT=ierr(3)) 171 174 ! 172 175 thd_ice_alloc = MAXVAL( ierr )
Note: See TracChangeset
for help on using the changeset viewer.
