- Timestamp:
- 2014-11-18T18:03:00+01:00 (9 years ago)
- Location:
- trunk/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 14 edited
-
ice.F90 (modified) (6 diffs)
-
limdyn.F90 (modified) (1 diff)
-
limitd_th.F90 (modified) (7 diffs)
-
limsbc.F90 (modified) (3 diffs)
-
limthd.F90 (modified) (5 diffs)
-
limthd_dh.F90 (modified) (33 diffs)
-
limthd_dif.F90 (modified) (45 diffs)
-
limthd_ent.F90 (modified) (1 diff)
-
limthd_lac.F90 (modified) (14 diffs)
-
limthd_sal.F90 (modified) (5 diffs)
-
limupdate1.F90 (modified) (3 diffs)
-
limupdate2.F90 (modified) (3 diffs)
-
limvar.F90 (modified) (5 diffs)
-
thd_ice.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r4869 r4872 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 | … … 346 346 !! * Old values of global variables 347 347 !!-------------------------------------------------------------------------- 348 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: old_v_s, old_v_i!: snow and ice volumes349 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: old_a_i, old_smv_i, old_oa_i !: ???350 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: old_e_s!: snow heat content351 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: old_e_i!: ice temperatures352 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 353 353 354 354 … … 487 487 ! * Old values of global variables 488 488 ii = ii + 1 489 ALLOCATE( old_v_s (jpi,jpj,jpl) , old_v_i (jpi,jpj,jpl) , old_e_s(jpi,jpj,nlay_s,jpl) , &490 & old_a_i (jpi,jpj,jpl) , old_smv_i(jpi,jpj,jpl) , old_e_i(jpi,jpj,jkmax ,jpl) , &491 & o ld_oa_i(jpi,jpj,jpl) , &492 & 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,jkmax ,jpl) , & 491 & oa_i_b (jpi,jpj,jpl) , & 492 & u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , STAT=ierr(ii) ) 493 493 494 494 ! * Increment of global variables -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r4863 r4872 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) -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r4869 r4872 166 166 REAL(wp), POINTER, DIMENSION(:,:,:) :: hL ! left boundary for the ITD for each thickness 167 167 REAL(wp), POINTER, DIMENSION(:,:,:) :: hR ! left boundary for the ITD for each thickness 168 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_ o! old ice thickness168 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_b ! old ice thickness 169 169 REAL(wp), POINTER, DIMENSION(:,:,:) :: dummy_es 170 170 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! local increment of ice area and volume … … 184 184 CALL wrk_alloc( jpi,jpj, zremap_flag ) ! integer 185 185 CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) ! integer 186 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 ) 187 187 CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice ) 188 188 CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) … … 221 221 zindb = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 222 222 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb 223 zindb = 1.0 - MAX( 0.0, SIGN( 1.0, - old_a_i(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes224 zht_i_ o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX( old_a_i(ji,jj,jl), epsi10 ) * zindb225 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) 226 226 END DO 227 227 END DO … … 268 268 ! 269 269 zhbnew(ii,ij,jl) = hi_max(jl) 270 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 271 271 !interpolate between adjacent category growth rates 272 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_ o(ii,ij,jl+1) - zht_i_o(ii,ij,jl) )273 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_ o(ii,ij,jl) )274 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 275 275 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 276 ELSEIF ( old_a_i(ii,ij,jl+1) > epsi10) THEN276 ELSEIF ( a_i_b(ii,ij,jl+1) > epsi10) THEN 277 277 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 278 278 ENDIF … … 337 337 !----------------------------------------------------------------------------------------------- 338 338 !- 7.1 g(h) for category 1 at start of time step 339 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_ o(:,:,klbnd), &339 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), & 340 340 & g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 341 341 & hR(:,:,klbnd), zremap_flag ) … … 362 362 ! Constrain new thickness <= ht_i 363 363 zdamax = a_i(ii,ij,klbnd) * & 364 (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 365 365 !ice area lost due to melting of thin ice 366 366 zda0 = MIN(zda0, zdamax) … … 484 484 CALL wrk_dealloc( jpi,jpj, zremap_flag ) ! integer 485 485 CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) ! integer 486 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 ) 487 487 CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice ) 488 488 CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r4787 r4872 117 117 CALL iom_put( "qsr_oce" , qsr(:,:) * pfrld(:,:) ) ! solar flux at ocean surface 118 118 CALL iom_put( "qns_oce" , qns(:,:) * pfrld(:,:) ) ! non-solar flux at ocean surface 119 CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) ) ! solar flux at ice surface120 CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) ) ! non-solar flux at ice surface121 CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) ) ! solar flux transmitted thru ice119 CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux at ice surface 120 CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! non-solar flux at ice surface 121 CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux transmitted thru ice 122 122 CALL iom_put( "qt_oce" , ( qsr(:,:) + qns(:,:) ) * pfrld(:,:) ) 123 CALL iom_put( "qt_ice" , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) * old_a_i(:,:,:), dim=3 ) )123 CALL iom_put( "qt_ice" , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) * a_i_b(:,:,:), dim=3 ) ) 124 124 125 125 ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) … … 139 139 !!!zfcm1 = qsr_tot(ji,jj) + ftr_ice(ji,jj) * ( 1._wp - pfrld(ji,jj) ) / ( 1._wp - zinda + zinda * iatte(ji,jj) ) 140 140 DO jl = 1, jpl 141 zfcm1 = zfcm1 - ( qsr_ice(ji,jj,jl) - ftr_ice(ji,jj,jl) ) * old_a_i(ji,jj,jl)141 zfcm1 = zfcm1 - ( qsr_ice(ji,jj,jl) - ftr_ice(ji,jj,jl) ) * a_i_b(ji,jj,jl) 142 142 END DO 143 143 ELSE … … 146 146 zfcm1 = pfrld(ji,jj) * qsr(ji,jj) 147 147 DO jl = 1, jpl 148 zfcm1 = zfcm1 + old_a_i(ji,jj,jl) * ftr_ice(ji,jj,jl)148 zfcm1 = zfcm1 + a_i_b(ji,jj,jl) * ftr_ice(ji,jj,jl) 149 149 END DO 150 150 ENDIF -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r4765 r4872 259 259 !------------------------- 260 260 261 CALL tab_2d_1d( nbpb, at_i_ b(1:nbpb), at_i , jpi, jpj, npb(1:nbpb) )262 CALL tab_2d_1d( nbpb, a_i_ b(1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )263 CALL tab_2d_1d( nbpb, ht_i_ b(1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )264 CALL tab_2d_1d( nbpb, ht_s_ b(1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) )265 266 CALL tab_2d_1d( nbpb, t_su_ b(1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) )267 CALL tab_2d_1d( nbpb, sm_i_ b(1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )261 CALL tab_2d_1d( nbpb, at_i_1d (1:nbpb), at_i , jpi, jpj, npb(1:nbpb) ) 262 CALL tab_2d_1d( nbpb, a_i_1d (1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 263 CALL tab_2d_1d( nbpb, ht_i_1d (1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 264 CALL tab_2d_1d( nbpb, ht_s_1d (1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 265 266 CALL tab_2d_1d( nbpb, t_su_1d (1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 267 CALL tab_2d_1d( nbpb, sm_i_1d (1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 268 268 DO jk = 1, nlay_s 269 CALL tab_2d_1d( nbpb, t_s_ b(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )270 CALL tab_2d_1d( nbpb, q_s_ b(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )269 CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 270 CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 271 271 END DO 272 272 DO jk = 1, nlay_i 273 CALL tab_2d_1d( nbpb, t_i_ b(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )274 CALL tab_2d_1d( nbpb, q_i_ b(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )275 CALL tab_2d_1d( nbpb, s_i_ b(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )273 CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 274 CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 275 CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 276 276 END DO 277 277 … … 287 287 ENDIF 288 288 CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 289 CALL tab_2d_1d( nbpb, t_bo_ b(1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) )289 CALL tab_2d_1d( nbpb, t_bo_1d (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) ) 290 290 CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip , jpi, jpj, npb(1:nbpb) ) 291 291 CALL tab_2d_1d( nbpb, fhtur_1d (1:nbpb), fhtur , jpi, jpj, npb(1:nbpb) ) … … 341 341 342 342 ! --- Ice enthalpy remapping --- ! 343 CALL lim_thd_ent( 1, nbpb, q_i_ b(1:nbpb,:) )343 CALL lim_thd_ent( 1, nbpb, q_i_1d(1:nbpb,:) ) 344 344 345 345 !---------------------------------! … … 357 357 !-------------------------------- 358 358 359 CALL tab_1d_2d( nbpb, at_i , npb, at_i_ b(1:nbpb) , jpi, jpj )360 CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_ b(1:nbpb) , jpi, jpj )361 CALL tab_1d_2d( nbpb, ht_s(:,:,jl) , npb, ht_s_ b(1:nbpb) , jpi, jpj )362 CALL tab_1d_2d( nbpb, a_i (:,:,jl) , npb, a_i_ b(1:nbpb) , jpi, jpj )363 CALL tab_1d_2d( nbpb, t_su(:,:,jl) , npb, t_su_ b(1:nbpb) , jpi, jpj )364 CALL tab_1d_2d( nbpb, sm_i(:,:,jl) , npb, sm_i_ b(1:nbpb) , jpi, jpj )359 CALL tab_1d_2d( nbpb, at_i , npb, at_i_1d (1:nbpb) , jpi, jpj ) 360 CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_1d (1:nbpb) , jpi, jpj ) 361 CALL tab_1d_2d( nbpb, ht_s(:,:,jl) , npb, ht_s_1d (1:nbpb) , jpi, jpj ) 362 CALL tab_1d_2d( nbpb, a_i (:,:,jl) , npb, a_i_1d (1:nbpb) , jpi, jpj ) 363 CALL tab_1d_2d( nbpb, t_su(:,:,jl) , npb, t_su_1d (1:nbpb) , jpi, jpj ) 364 CALL tab_1d_2d( nbpb, sm_i(:,:,jl) , npb, sm_i_1d (1:nbpb) , jpi, jpj ) 365 365 DO jk = 1, nlay_s 366 CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_ b(1:nbpb,jk), jpi, jpj)367 CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_ b(1:nbpb,jk), jpi, jpj)366 CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d (1:nbpb,jk), jpi, jpj) 367 CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d (1:nbpb,jk), jpi, jpj) 368 368 END DO 369 369 DO jk = 1, nlay_i 370 CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_ b(1:nbpb,jk), jpi, jpj)371 CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_ b(1:nbpb,jk), jpi, jpj)372 CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_ b(1:nbpb,jk), jpi, jpj)370 CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d (1:nbpb,jk), jpi, jpj) 371 CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d (1:nbpb,jk), jpi, jpj) 372 CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d (1:nbpb,jk), jpi, jpj) 373 373 END DO 374 374 CALL tab_1d_2d( nbpb, qlead , npb, qlead_1d (1:nbpb) , jpi, jpj ) … … 507 507 DO jk = 1, nlay_i 508 508 DO ji = kideb, kiut 509 ztmelts = -tmut * s_i_ b(ji,jk) + rtt509 ztmelts = -tmut * s_i_1d(ji,jk) + rtt 510 510 ! Conversion q(S,T) -> T (second order equation) 511 511 zaaa = cpic 512 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_ b(ji,jk) / rhoic - lfus512 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_1d(ji,jk) / rhoic - lfus 513 513 zccc = lfus * ( ztmelts - rtt ) 514 514 zdiscrim = SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) 515 t_i_ b(ji,jk) = rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa )515 t_i_1d(ji,jk) = rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) 516 516 517 517 ! mask temperature 518 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_ b(ji) ) )519 t_i_ b(ji,jk) = zswitch * t_i_b(ji,jk) + ( 1._wp - zswitch ) * rtt518 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 519 t_i_1d(ji,jk) = zswitch * t_i_1d(ji,jk) + ( 1._wp - zswitch ) * rtt 520 520 END DO 521 521 END DO -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r4688 r4872 156 156 DO jk = 1, nlay_i 157 157 DO ji = kideb, kiut 158 h_i_old (ji,jk) = ht_i_ b(ji) / REAL( nlay_i )159 qh_i_old(ji,jk) = q_i_ b(ji,jk) * h_i_old(ji,jk)158 h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 159 qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk) 160 160 ENDDO 161 161 ENDDO … … 166 166 ! 167 167 DO ji = kideb, kiut 168 zinda = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_ b(ji) ) )168 zinda = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 169 169 ztmelts = zinda * rtt + ( 1._wp - zinda ) * rtt 170 170 … … 172 172 zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 173 173 174 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_ b(ji) - ztmelts ) )174 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts ) ) 175 175 zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 176 176 END DO … … 182 182 !------------------------------------------------------------------------------! 183 183 DO ji = kideb, kiut 184 IF( t_s_ b(ji,1) > rtt ) THEN !!! Internal melting184 IF( t_s_1d(ji,1) > rtt ) THEN !!! Internal melting 185 185 ! Contribution to heat flux to the ocean [W.m-2], < 0 186 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_ b(ji,1) * ht_s_b(ji) * a_i_b(ji) * r1_rdtice186 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 187 187 ! Contribution to mass flux 188 wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_ b(ji) * a_i_b(ji) * r1_rdtice188 wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 189 189 ! updates 190 ht_s_ b(ji) = 0._wp191 q_s_ b(ji,1) = 0._wp192 t_s_ b(ji,1) = rtt190 ht_s_1d(ji) = 0._wp 191 q_s_1d (ji,1) = 0._wp 192 t_s_1d (ji,1) = rtt 193 193 END IF 194 194 END DO … … 199 199 ! 200 200 DO ji = kideb, kiut 201 zh_s(ji) = ht_s_ b(ji) / REAL( nlay_s )201 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 202 202 END DO 203 203 ! 204 204 DO jk = 1, nlay_s 205 205 DO ji = kideb, kiut 206 zqh_s(ji) = zqh_s(ji) + q_s_ b(ji,jk) * zh_s(ji)206 zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji) 207 207 END DO 208 208 END DO … … 210 210 DO jk = 1, nlay_i 211 211 DO ji = kideb, kiut 212 zh_i(ji,jk) = ht_i_ b(ji) / REAL( nlay_i )213 zqh_i(ji) = zqh_i(ji) + q_i_ b(ji,jk) * zh_i(ji,jk)212 zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 213 zqh_i(ji) = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk) 214 214 END DO 215 215 END DO … … 238 238 !----------- 239 239 ! thickness change 240 zcoeff = ( 1._wp - ( 1._wp - at_i_ b(ji) )**betas ) / at_i_b(ji)240 zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**betas ) / at_i_1d(ji) 241 241 zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn 242 242 ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) … … 244 244 IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 245 245 ! heat flux from snow precip (>0, W.m-2) 246 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_ b(ji) * zqprec(ji) * r1_rdtice246 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 247 247 ! mass flux, <0 248 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_ b(ji) * zdh_s_pre(ji) * r1_rdtice248 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice 249 249 ! update thickness 250 ht_s_ b (ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_pre(ji) )250 ht_s_1d (ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) 251 251 252 252 !--------------------- … … 259 259 zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting 260 260 ! heat used to melt snow (W.m-2, >0) 261 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_ b(ji) * zqprec(ji) * r1_rdtice261 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 262 262 ! snow melting only = water into the ocean (then without snow precip), >0 263 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_ b(ji) * zdh_s_mel(ji) * r1_rdtice263 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_mel(ji) * r1_rdtice 264 264 265 265 ! updates available heat + thickness 266 266 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) ) 267 ht_s_ b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_mel(ji) )268 zh_s (ji) = ht_s_ b(ji) / REAL( nlay_s )267 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) ) 268 zh_s (ji) = ht_s_1d(ji) / REAL( nlay_s ) 269 269 270 270 ENDIF … … 276 276 DO ji = kideb, kiut 277 277 ! thickness change 278 zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_ b(ji) ) )279 zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_ b(ji,jk) + epsi20 ) )280 zdeltah (ji,jk) = - zindh * zindq * zq_su(ji) / MAX( q_s_ b(ji,jk), epsi20 )278 zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 279 zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_1d(ji,jk) + epsi20 ) ) 280 zdeltah (ji,jk) = - zindh * zindq * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 ) 281 281 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting 282 282 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 283 283 ! heat used to melt snow(W.m-2, >0) 284 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_ b(ji) * q_s_b(ji,jk) * r1_rdtice284 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice 285 285 ! snow melting only = water into the ocean (then without snow precip) 286 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice286 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 287 287 288 288 ! updates available heat + thickness 289 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_ b(ji,jk) )290 ht_s_ b(ji) = MAX( 0._wp , ht_s_b(ji) + zdeltah(ji,jk) )289 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 290 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 291 291 292 292 END DO … … 305 305 ! forced mode: snow thickness change due to sublimation 306 306 DO ji = kideb, kiut 307 zdh_s_sub(ji) = MAX( - ht_s_ b(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice )307 zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) 308 308 ! Heat flux by sublimation [W.m-2], < 0 309 309 ! sublimate first snow that had fallen, then pre-existing snow 310 310 zcoeff = ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) * zqprec(ji) + & 311 & ( 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) ) &312 & * a_i_ b(ji) * r1_rdtice311 & ( 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) ) & 312 & * a_i_1d(ji) * r1_rdtice 313 313 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff 314 314 ! Mass flux by sublimation 315 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_ b(ji) * zdh_s_sub(ji) * r1_rdtice315 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 316 316 ! new snow thickness 317 ht_s_ b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_sub(ji) )317 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 318 318 END DO 319 319 ENDIF … … 322 322 DO ji = kideb, kiut 323 323 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 324 zh_s(ji) = ht_s_ b(ji) / REAL( nlay_s )324 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 325 325 END DO ! ji 326 326 … … 332 332 DO jk = 1, nlay_s 333 333 DO ji = kideb,kiut 334 zindh = MAX( 0._wp , SIGN( 1._wp, - ht_s_ b(ji) + epsi20 ) )335 q_s_ b(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_b(ji), epsi20 ) * &334 zindh = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) + epsi20 ) ) 335 q_s_1d(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_1d(ji), epsi20 ) * & 336 336 & ( ( MAX( 0._wp, dh_s_tot(ji) ) ) * zqprec(ji) + & 337 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_ b(ji) ) * rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) )338 zq_s(ji) = zq_s(ji) + q_s_ b(ji,jk)337 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) ) 338 zq_s(ji) = zq_s(ji) + q_s_1d(ji,jk) 339 339 END DO 340 340 END DO … … 346 346 DO jk = 1, nlay_i 347 347 DO ji = kideb, kiut 348 zEi = - q_i_ b(ji,jk) / rhoic ! Specific enthalpy of layer k [J/kg, <0]349 350 ztmelts = - tmut * s_i_ b(ji,jk) + rtt ! Melting point of layer k [K]348 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of layer k [J/kg, <0] 349 350 ztmelts = - tmut * s_i_1d(ji,jk) + rtt ! Melting point of layer k [K] 351 351 352 352 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] … … 368 368 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 369 369 370 ! Contribution to salt flux (clem: using sm_i_ b and not s_i_b(jk) is ok)371 sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_ b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice370 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 371 sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 372 372 373 373 ! Contribution to heat flux [W.m-2], < 0 374 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice374 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 375 375 376 376 ! Total heat flux used in this process [W.m-2], > 0 377 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_ b(ji) * zdE * r1_rdtice377 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 378 378 379 379 ! Contribution to mass flux 380 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice380 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 381 381 382 382 ! record which layers have disappeared (for bottom melting) … … 388 388 389 389 ! update heat content (J.m-2) and layer thickness 390 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_ b(ji,jk)390 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 391 391 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 392 392 END DO … … 394 394 ! update ice thickness 395 395 DO ji = kideb, kiut 396 ht_i_ b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_surf(ji) )396 ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) ) 397 397 END DO 398 398 … … 424 424 !clem debug. Just to be sure that enthalpy at nlay_i+1 is null 425 425 DO ji = kideb, kiut 426 q_i_ b(ji,nlay_i+1) = 0._wp426 q_i_1d(ji,nlay_i+1) = 0._wp 427 427 END DO 428 428 … … 446 446 447 447 s_i_new(ji) = zswitch_sal * zfracs * sss_m(ii,ij) & ! New ice salinity 448 + ( 1. - zswitch_sal ) * sm_i_ b(ji)448 + ( 1. - zswitch_sal ) * sm_i_1d(ji) 449 449 ! New ice growth 450 450 ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K) 451 451 452 zt_i_new = zswitch_sal * t_bo_ b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)452 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 453 453 454 454 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) … … 456 456 & + rcp * ( ztmelts-rtt ) 457 457 458 zEw = rcp * ( t_bo_ b(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)458 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 459 459 460 460 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) … … 462 462 dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 463 463 464 q_i_ b(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0)464 q_i_1d(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0) 465 465 466 466 ENDIF ! fc_bo_i … … 477 477 ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K) 478 478 479 zt_i_new = zswitch_sal * t_bo_ b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)479 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 480 480 481 481 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) … … 483 483 & + rcp * ( ztmelts-rtt ) 484 484 485 zEw = rcp * ( t_bo_ b(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)485 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 486 486 487 487 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 488 488 489 489 ! Contribution to heat flux to the ocean [W.m-2], >0 490 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice490 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 491 491 492 492 ! Total heat flux used in this process [W.m-2], <0 493 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_ b(ji) * zdE * r1_rdtice493 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 494 494 495 495 ! Contribution to salt flux, <0 496 sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_ b(ji) * zfmdt * r1_rdtice496 sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_1d(ji) * zfmdt * r1_rdtice 497 497 498 498 ! Contribution to mass flux, <0 499 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_ b(ji) * dh_i_bott(ji) * r1_rdtice499 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice 500 500 501 501 ! update heat content (J.m-2) and layer thickness 502 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)502 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) 503 503 h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji) 504 504 ENDIF … … 513 513 IF( zf_tt(ji) >= 0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared from surface melting 514 514 515 ztmelts = - tmut * s_i_ b(ji,jk) + rtt ! Melting point of layer jk (K)516 517 IF( t_i_ b(ji,jk) >= ztmelts ) THEN !!! Internal melting515 ztmelts = - tmut * s_i_1d(ji,jk) + rtt ! Melting point of layer jk (K) 516 517 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 518 518 zintermelt(ji) = 1._wp 519 519 520 zEi = - q_i_ b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)521 522 !!zEw = rcp * ( t_i_ b(ji,jk) - rtt ) ! Specific enthalpy of meltwater at T = t_i_b(J/kg, <0)520 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 521 522 !!zEw = rcp * ( t_i_1d(ji,jk) - rtt ) ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) 523 523 524 524 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) … … 533 533 534 534 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 535 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_ b(ji) * zEi * r1_rdtice536 537 ! Contribution to salt flux (clem: using sm_i_ b and not s_i_b(jk) is ok)538 sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_ b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice535 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 536 537 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 538 sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 539 539 540 540 ! Contribution to mass flux 541 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice541 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 542 542 543 543 ! update heat content (J.m-2) and layer thickness 544 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_ b(ji,jk)544 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 545 545 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 546 546 547 547 ELSE !!! Basal melting 548 548 549 zEi = - q_i_ b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)549 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 550 550 551 551 zEw = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0) … … 568 568 569 569 ! Contribution to heat flux to the ocean [W.m-2], <0 570 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice571 572 ! Contribution to salt flux (clem: using sm_i_ b and not s_i_b(jk) is ok)573 sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_ b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice570 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 571 572 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 573 sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 574 574 575 575 ! Total heat flux used in this process [W.m-2], >0 576 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_ b(ji) * zdE * r1_rdtice576 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 577 577 578 578 ! Contribution to mass flux 579 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice579 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 580 580 581 581 ! update heat content (J.m-2) and layer thickness 582 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_ b(ji,jk)582 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 583 583 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 584 584 ENDIF … … 603 603 ! 604 604 ! ! excessive energy is sent to lateral ablation 605 ! zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_ b(ji) - epsi20 ) )606 ! zq_1cat(ji) = zinda * rhoic * lfus * at_i_ b(ji) / MAX( 1._wp - at_i_b(ji) , epsi20 ) * zdvres ! J.m-2 >=0605 ! zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_1d(ji) - epsi20 ) ) 606 ! zq_1cat(ji) = zinda * rhoic * lfus * at_i_1d(ji) / MAX( 1._wp - at_i_1d(ji) , epsi20 ) * zdvres ! J.m-2 >=0 607 607 ! 608 608 ! ! correct salt and mass fluxes 609 ! 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 approximation610 ! wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_ b(ji) * zdvres * r1_rdtice609 ! 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 610 ! wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdvres * r1_rdtice 611 611 ! ENDIF 612 612 ! END DO … … 617 617 !------------------------------------------- 618 618 DO ji = kideb, kiut 619 ht_i_ b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_bott(ji) )619 ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) ) 620 620 END DO 621 621 … … 628 628 DO ji = kideb, kiut 629 629 zq_rema(ji) = zq_su(ji) + zq_bo(ji) 630 ! zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_ b(ji) ) ) ! =1 if snow630 ! zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) ! =1 if snow 631 631 ! zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zq_s(ji) + epsi20 ) ) 632 632 ! zdeltah (ji,1) = - zindh * zindq * zq_rema(ji) / MAX( zq_s(ji), epsi20 ) 633 ! zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_ b(ji) ) ) ! bound melting633 ! zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 634 634 ! zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1) 635 635 ! dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1) 636 ! ht_s_ b (ji) = ht_s_b(ji) + zdeltah(ji,1)636 ! ht_s_1d (ji) = ht_s_1d(ji) + zdeltah(ji,1) 637 637 ! 638 638 ! zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji) ! update available heat (J.m-2) 639 639 ! ! heat used to melt snow 640 ! hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_ b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0)640 ! hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0) 641 641 ! ! Contribution to mass flux 642 ! wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_ b(ji) * zdeltah(ji,1) * r1_rdtice642 ! wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 643 643 ! 644 644 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 645 645 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 646 hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_ b(ji) ) * r1_rdtice646 hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 647 647 648 648 IF( ln_nicep .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) … … 657 657 DO ji = kideb, kiut 658 658 ! 659 dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_ b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic ) )660 661 ht_i_ b(ji) = ht_i_b(ji) + dh_snowice(ji)662 ht_s_ b(ji) = ht_s_b(ji) - dh_snowice(ji)659 dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic ) ) 660 661 ht_i_1d(ji) = ht_i_1d(ji) + dh_snowice(ji) 662 ht_s_1d(ji) = ht_s_1d(ji) - dh_snowice(ji) 663 663 664 664 ! Salinity of snow ice 665 665 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 666 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_ b(ji)666 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 667 667 668 668 ! entrapment during snow ice formation 669 669 ! new salinity difference stored (to be used in limthd_ent.F90) 670 670 IF ( num_sal == 2 ) THEN 671 zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_ b(ji) - epsi10 ) )671 zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) ) 672 672 ! salinity dif due to snow-ice formation 673 dsm_i_si_1d(ji) = ( zs_snic - sm_i_ b(ji) ) * dh_snowice(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch673 dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi10 ) * zswitch 674 674 ! salinity dif due to bottom growth 675 675 IF ( zf_tt(ji) < 0._wp ) THEN 676 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_ b(ji) ) * dh_i_bott(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch676 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi10 ) * zswitch 677 677 ENDIF 678 678 ENDIF … … 686 686 687 687 ! Contribution to heat flux 688 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice688 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 689 689 690 690 ! Contribution to salt flux 691 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_ b(ji) * zfmdt * r1_rdtice691 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice 692 692 693 693 ! Contribution to mass flux 694 694 ! All snow is thrown in the ocean, and seawater is taken to replace the volume 695 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_ b(ji) * dh_snowice(ji) * rhoic * r1_rdtice696 wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_ b(ji) * dh_snowice(ji) * rhosn * r1_rdtice695 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice 696 wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice 697 697 698 698 ! update heat content (J.m-2) and layer thickness 699 qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_ b(ji,1) + zfmdt * zEw699 qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw 700 700 h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 701 701 702 702 ! Total ablation (to debug) 703 IF( ht_i_ b(ji) <= 0._wp ) a_i_b(ji) = 0._wp703 IF( ht_i_1d(ji) <= 0._wp ) a_i_1d(ji) = 0._wp 704 704 705 705 END DO !ji … … 711 711 !clem bug: we should take snow into account here 712 712 DO ji = kideb, kiut 713 zindh = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_ b(ji) ) )714 t_su_ b(ji) = zindh * t_su_b(ji) + ( 1.0 - zindh ) * rtt713 zindh = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 714 t_su_1d(ji) = zindh * t_su_1d(ji) + ( 1.0 - zindh ) * rtt 715 715 END DO ! ji 716 716 … … 718 718 DO ji = kideb,kiut 719 719 ! mask enthalpy 720 zinda = MAX( 0._wp , SIGN( 1._wp, - ht_s_ b(ji) ) )721 q_s_ b(ji,jk) = ( 1.0 - zinda ) * q_s_b(ji,jk)722 ! recalculate t_s_ b from q_s_b723 t_s_ b(ji,jk) = rtt + ( 1._wp - zinda ) * ( - q_s_b(ji,jk) / ( rhosn * cpic ) + lfus / cpic )720 zinda = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) ) ) 721 q_s_1d(ji,jk) = ( 1.0 - zinda ) * q_s_1d(ji,jk) 722 ! recalculate t_s_1d from q_s_1d 723 t_s_1d(ji,jk) = rtt + ( 1._wp - zinda ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 724 724 END DO 725 725 END DO -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r4870 r4872 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 : … … 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 … … 151 151 ! 152 152 CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow ) 153 CALL wrk_alloc( jpij, ztfs, ztsu old, ztsuoldit, zh_i, zh_s, zfsw )153 CALL wrk_alloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 154 154 CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 155 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)156 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, zts old, zeta_s, ztstemp, z_s, kjstart=0)155 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) 156 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0) 157 157 CALL wrk_alloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis ) 158 158 CALL wrk_alloc( jpij, jkmax+2, 3, ztrid ) … … 163 163 zdq(:) = 0._wp ; zq_ini(:) = 0._wp 164 164 DO ji = kideb, kiut 165 zq_ini(ji) = ( SUM( q_i_ b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + &166 & SUM( q_s_ b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )165 zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) + & 166 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 167 167 END DO 168 168 … … 173 173 DO ji = kideb , kiut 174 174 ! is there snow or not 175 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_ b(ji) ) ) )175 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) ) ) 176 176 ! surface temperature of fusion 177 177 ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 178 178 ! layer thickness 179 zh_i(ji) = ht_i_ b(ji) / REAL( nlay_i )180 zh_s(ji) = ht_s_ b(ji) / REAL( nlay_s )179 zh_i(ji) = ht_i_1d(ji) / REAL( nlay_i ) 180 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 181 181 END DO 182 182 … … 190 190 DO jk = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 191 191 DO ji = kideb , kiut 192 z_s(ji,jk) = z_s(ji,jk-1) + ht_s_ b(ji) / REAL( nlay_s )192 z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) / REAL( nlay_s ) 193 193 END DO 194 194 END DO … … 196 196 DO jk = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 197 197 DO ji = kideb , kiut 198 z_i(ji,jk) = z_i(ji,jk-1) + ht_i_ b(ji) / REAL( nlay_i )198 z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) / REAL( nlay_i ) 199 199 END DO 200 200 END DO … … 217 217 DO ji = kideb , kiut 218 218 ! switches 219 isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_ b(ji) ) ) )219 isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) ) 220 220 ! hs > 0, isnow = 1 221 221 zhsu (ji) = hnzst ! threshold for the computation of i0 222 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_ b(ji) / zhsu(ji) ) )222 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu(ji) ) ) 223 223 224 224 i0(ji) = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) … … 227 227 ! a function of the cloud cover 228 228 ! 229 !i0(ji) = (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_ b(ji)+10.0)229 !i0(ji) = (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_1d(ji)+10.0) 230 230 !formula used in Cice 231 231 END DO … … 272 272 273 273 DO ji = kideb, kiut ! Radiation transmitted below the ice 274 !!!ftr_ice_1d(ji) = ftr_ice_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_ b(ji) / at_i_b(ji) ! clem modif274 !!!ftr_ice_1d(ji) = ftr_ice_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_1d(ji) / at_i_1d(ji) ! clem modif 275 275 ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 276 276 END DO … … 282 282 ! 283 283 DO ji = kideb, kiut ! Old surface temperature 284 ztsu old (ji) = t_su_b(ji) ! temperature at the beg of iter pr.285 ztsu oldit(ji) = t_su_b(ji) ! temperature at the previous iter286 t_su_ b (ji) = MIN( t_su_b(ji), ztfs(ji) - ztsu_err ) ! necessary284 ztsub (ji) = t_su_1d(ji) ! temperature at the beg of iter pr. 285 ztsubit(ji) = t_su_1d(ji) ! temperature at the previous iter 286 t_su_1d (ji) = MIN( t_su_1d(ji), ztfs(ji) - ztsu_err ) ! necessary 287 287 zerrit (ji) = 1000._wp ! initial value of error 288 288 END DO … … 290 290 DO jk = 1, nlay_s ! Old snow temperature 291 291 DO ji = kideb , kiut 292 zts old(ji,jk) = t_s_b(ji,jk)292 ztsb(ji,jk) = t_s_1d(ji,jk) 293 293 END DO 294 294 END DO … … 296 296 DO jk = 1, nlay_i ! Old ice temperature 297 297 DO ji = kideb , kiut 298 zti old(ji,jk) = t_i_b(ji,jk)298 ztib(ji,jk) = t_i_1d(ji,jk) 299 299 END DO 300 300 END DO … … 313 313 IF( thcon_i_swi == 0 ) THEN ! Untersteiner (1964) formula 314 314 DO ji = kideb , kiut 315 ztcond_i(ji,0) = rcdic + zbeta*s_i_ b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt)315 ztcond_i(ji,0) = rcdic + zbeta*s_i_1d(ji,1) / MIN(-epsi10,t_i_1d(ji,1)-rtt) 316 316 ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin) 317 317 END DO 318 318 DO jk = 1, nlay_i-1 319 319 DO ji = kideb , kiut 320 ztcond_i(ji,jk) = rcdic + zbeta*( s_i_ b(ji,jk) + s_i_b(ji,jk+1) ) / &321 MIN(-2.0_wp * epsi10, t_i_ b(ji,jk)+t_i_b(ji,jk+1) - 2.0_wp * rtt)320 ztcond_i(ji,jk) = rcdic + zbeta*( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / & 321 MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt) 322 322 ztcond_i(ji,jk) = MAX(ztcond_i(ji,jk),zkimin) 323 323 END DO … … 327 327 IF( thcon_i_swi == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 328 328 DO ji = kideb , kiut 329 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_ b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt ) &330 & - 0.011_wp * ( t_i_ b(ji,1) - rtt )329 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1)-rtt ) & 330 & - 0.011_wp * ( t_i_1d(ji,1) - rtt ) 331 331 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 332 332 END DO … … 334 334 DO ji = kideb , kiut 335 335 ztcond_i(ji,jk) = rcdic + & 336 & 0.090_wp * ( s_i_ b(ji,jk) + s_i_b(ji,jk+1) ) &337 & / MIN(-2.0_wp * epsi10, t_i_ b(ji,jk)+t_i_b(ji,jk+1) - 2.0_wp * rtt) &338 & - 0.0055_wp* ( t_i_ b(ji,jk) + t_i_b(ji,jk+1) - 2.0*rtt )336 & 0.090_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) & 337 & / MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt) & 338 & - 0.0055_wp* ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0*rtt ) 339 339 ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 340 340 END DO 341 341 END DO 342 342 DO ji = kideb , kiut 343 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_ b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt) &344 & - 0.011_wp * ( t_bo_ b(ji) - rtt )343 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN(-epsi10,t_bo_1d(ji)-rtt) & 344 & - 0.011_wp * ( t_bo_1d(ji) - rtt ) 345 345 ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 346 346 END DO … … 389 389 DO jk = 1, nlay_i 390 390 DO ji = kideb , kiut 391 ztitemp(ji,jk) = t_i_ b(ji,jk)392 zspeche_i(ji,jk) = cpic + zgamma*s_i_ b(ji,jk)/ &393 MAX((t_i_ b(ji,jk)-rtt)*(ztiold(ji,jk)-rtt),epsi10)391 ztitemp(ji,jk) = t_i_1d(ji,jk) 392 zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ & 393 MAX((t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt),epsi10) 394 394 zeta_i(ji,jk) = rdt_ice / MAX(rhoic*zspeche_i(ji,jk)*zh_i(ji), & 395 395 epsi10) … … 399 399 DO jk = 1, nlay_s 400 400 DO ji = kideb , kiut 401 ztstemp(ji,jk) = t_s_ b(ji,jk)401 ztstemp(ji,jk) = t_s_1d(ji,jk) 402 402 zeta_s(ji,jk) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10) 403 403 END DO … … 410 410 DO ji = kideb , kiut 411 411 ! update of the non solar flux according to the update in T_su 412 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_ b(ji) - ztsuoldit(ji) )412 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) ) 413 413 414 414 ! update incoming flux … … 448 448 zkappa_i(ji,jk)) 449 449 ztrid(ji,numeq,3) = - zeta_i(ji,jk)*zkappa_i(ji,jk) 450 zindterm(ji,numeq) = zti old(ji,jk) + zeta_i(ji,jk)* &450 zindterm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk)* & 451 451 zradab_i(ji,jk) 452 452 END DO … … 460 460 + zkappa_i(ji,nlay_i-1) ) 461 461 ztrid(ji,numeq,3) = 0.0 462 zindterm(ji,numeq) = zti old(ji,nlay_i) + zeta_i(ji,nlay_i)* &462 zindterm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 463 463 ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 464 * t_bo_ b(ji) )464 * t_bo_1d(ji) ) 465 465 ENDDO 466 466 467 467 468 468 DO ji = kideb , kiut 469 IF ( ht_s_ b(ji).gt.0.0 ) THEN469 IF ( ht_s_1d(ji).gt.0.0 ) THEN 470 470 ! 471 471 !------------------------------------------------------------------------------| … … 480 480 zkappa_s(ji,jk) ) 481 481 ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 482 zindterm(ji,numeq) = zts old(ji,jk) + zeta_s(ji,jk)* &482 zindterm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk)* & 483 483 zradab_s(ji,jk) 484 484 END DO … … 488 488 ztrid(ji,nlay_s+2,3) = 0.0 489 489 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 490 t_bo_ b(ji)490 t_bo_1d(ji) 491 491 ENDIF 492 492 493 IF ( t_su_ b(ji) .LT. rtt ) THEN493 IF ( t_su_1d(ji) .LT. rtt ) THEN 494 494 495 495 !------------------------------------------------------------------------------| … … 504 504 ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 505 505 ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 506 zindterm(ji,1) = dzf(ji)*t_su_ b(ji) - zf(ji)506 zindterm(ji,1) = dzf(ji)*t_su_1d(ji) - zf(ji) 507 507 508 508 !!first layer of snow equation … … 510 510 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s) 511 511 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 512 zindterm(ji,2) = zts old(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)512 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 513 513 514 514 ELSE … … 527 527 zkappa_s(ji,0) * zg1s ) 528 528 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 529 zindterm(ji,2) = zts old(ji,1) + zeta_s(ji,1) * &529 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * & 530 530 ( zradab_s(ji,1) + & 531 zkappa_s(ji,0) * zg1s * t_su_ b(ji) )531 zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 532 532 ENDIF 533 533 ELSE … … 537 537 !------------------------------------------------------------------------------| 538 538 ! 539 IF (t_su_ b(ji) .LT. rtt) THEN539 IF (t_su_1d(ji) .LT. rtt) THEN 540 540 ! 541 541 !------------------------------------------------------------------------------| … … 551 551 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1 552 552 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 553 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_ b(ji) - zf(ji)553 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji) 554 554 555 555 !!first layer of ice equation … … 558 558 + zkappa_i(ji,0) * zg1 ) 559 559 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1) 560 zindterm(ji,numeqmin(ji)+1)= zti old(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)560 zindterm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1) 561 561 562 562 !!case of only one layer in the ice (surface & ice equations are altered) … … 571 571 ztrid(ji,numeqmin(ji)+1,3) = 0.0 572 572 573 zindterm(ji,numeqmin(ji)+1) = zti old(ji,1) + zeta_i(ji,1)* &574 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_ b(ji) )573 zindterm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1)* & 574 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) ) 575 575 ENDIF 576 576 … … 591 591 zg1) 592 592 ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 593 zindterm(ji,numeqmin(ji)) = zti old(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &594 zkappa_i(ji,0) * zg1 * t_su_ b(ji) )593 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 594 zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 595 595 596 596 !!case of only one layer in the ice (surface & ice equations are altered) … … 600 600 zkappa_i(ji,1)) 601 601 ztrid(ji,numeqmin(ji),3) = 0.0 602 zindterm(ji,numeqmin(ji)) = zti old(ji,1) + zeta_i(ji,1)* &603 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_ b(ji)) &604 + t_su_ b(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0602 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)* & 603 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) & 604 + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0 605 605 ENDIF 606 606 … … 642 642 DO ji = kideb , kiut 643 643 ! ice temperatures 644 t_i_ b(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))644 t_i_1d(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 645 645 END DO 646 646 … … 648 648 DO ji = kideb , kiut 649 649 jk = numeq - nlay_s - 1 650 t_i_ b(ji,jk) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &651 t_i_ b(ji,jk+1))/zdiagbis(ji,numeq)650 t_i_1d(ji,jk) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* & 651 t_i_1d(ji,jk+1))/zdiagbis(ji,numeq) 652 652 END DO 653 653 END DO … … 655 655 DO ji = kideb , kiut 656 656 ! snow temperatures 657 IF (ht_s_ b(ji).GT.0._wp) &658 t_s_ b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &659 * t_i_ b(ji,1))/zdiagbis(ji,nlay_s+1) &660 * MAX(0.0,SIGN(1.0,ht_s_ b(ji)))657 IF (ht_s_1d(ji).GT.0._wp) & 658 t_s_1d(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 659 * t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) & 660 * MAX(0.0,SIGN(1.0,ht_s_1d(ji))) 661 661 662 662 ! surface temperature 663 isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_ b(ji) ) ) )664 ztsu oldit(ji) = t_su_b(ji)665 IF( t_su_ b(ji) < ztfs(ji) ) &666 t_su_ b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1) &667 & + REAL( 1 - isnow(ji) )*t_i_ b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))663 isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_1d(ji) ) ) ) 664 ztsubit(ji) = t_su_1d(ji) 665 IF( t_su_1d(ji) < ztfs(ji) ) & 666 t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1) & 667 & + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 668 668 END DO 669 669 ! … … 675 675 ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd 676 676 DO ji = kideb , kiut 677 t_su_ b(ji) = MAX( MIN( t_su_b(ji) , ztfs(ji) ) , 190._wp )678 zerrit(ji) = ABS( t_su_ b(ji) - ztsuoldit(ji) )677 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , ztfs(ji) ) , 190._wp ) 678 zerrit(ji) = ABS( t_su_1d(ji) - ztsubit(ji) ) 679 679 END DO 680 680 681 681 DO jk = 1, nlay_s 682 682 DO ji = kideb , kiut 683 t_s_ b(ji,jk) = MAX( MIN( t_s_b(ji,jk), rtt ), 190._wp )684 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_ b(ji,jk) - ztstemp(ji,jk)))683 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rtt ), 190._wp ) 684 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_1d(ji,jk) - ztstemp(ji,jk))) 685 685 END DO 686 686 END DO … … 688 688 DO jk = 1, nlay_i 689 689 DO ji = kideb , kiut 690 ztmelt_i = -tmut * s_i_ b(ji,jk) + rtt691 t_i_ b(ji,jk) = MAX(MIN(t_i_b(ji,jk),ztmelt_i), 190._wp)692 zerrit(ji) = MAX(zerrit(ji),ABS(t_i_ b(ji,jk) - ztitemp(ji,jk)))690 ztmelt_i = -tmut * s_i_1d(ji,jk) + rtt 691 t_i_1d(ji,jk) = MAX(MIN(t_i_1d(ji,jk),ztmelt_i), 190._wp) 692 zerrit(ji) = MAX(zerrit(ji),ABS(t_i_1d(ji,jk) - ztitemp(ji,jk))) 693 693 END DO 694 694 END DO … … 715 715 DO ji = kideb, kiut 716 716 ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux) 717 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) ) )717 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) ) ) 718 718 ! ! surface ice conduction flux 719 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_ b(ji) ) ) )720 fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_ b(ji,1) - t_su_b(ji)) &721 & - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_ b(ji,1) - t_su_b(ji))719 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) ) 720 fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & 721 & - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji)) 722 722 ! ! bottom ice conduction flux 723 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_ b(ji) - t_i_b(ji,nlay_i)) )723 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 724 724 END DO 725 725 … … 728 728 !----------------------------------------- 729 729 DO ji = kideb, kiut 730 IF( t_su_ b(ji) < rtt ) THEN ! case T_su < 0degC730 IF( t_su_1d(ji) < rtt ) THEN ! case T_su < 0degC 731 731 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 732 & ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_ b(ji)732 & ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 733 733 ELSE ! case T_su = 0degC 734 734 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 735 & ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_ b(ji)735 & ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 736 736 ENDIF 737 737 END DO … … 742 742 ! --- diag error on heat diffusion - PART 2 --- ! 743 743 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 ) )744 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) + & 745 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 746 746 zhfx_err = ( 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 * a_i_ b(ji)747 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err * a_i_1d(ji) 748 748 ! --- correction of qns_ice and surface conduction flux --- ! 749 749 qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err … … 751 751 ! --- Heat flux at the ice surface in W.m-2 --- ! 752 752 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 753 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_ b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) )753 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_1d(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 754 754 END DO 755 755 756 756 ! 757 757 CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 758 CALL wrk_dealloc( jpij, ztfs, ztsu old, ztsuoldit, zh_i, zh_s, zfsw )758 CALL wrk_dealloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 759 759 CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 760 760 CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, & 761 & zti old, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 )762 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, zts old, zeta_s, ztstemp, z_s, kjstart = 0 )761 & ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 762 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 763 763 CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis ) 764 764 CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid ) … … 783 783 DO jk = 1, nlay_i ! Sea ice energy of melting 784 784 DO ji = kideb, kiut 785 ztmelts = - tmut * s_i_ b(ji,jk) + rtt786 zindb = MAX( 0._wp , SIGN( 1._wp , -(t_i_ b(ji,jk) - rtt) - epsi10 ) )787 q_i_ b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) ) &788 & + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_ b(ji,jk)-rtt, -epsi10 ) ) &785 ztmelts = - tmut * s_i_1d(ji,jk) + rtt 786 zindb = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi10 ) ) 787 q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) & 788 & + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_1d(ji,jk)-rtt, -epsi10 ) ) & 789 789 & - rcp * ( ztmelts-rtt ) ) 790 790 END DO … … 792 792 DO jk = 1, nlay_s ! Snow energy of melting 793 793 DO ji = kideb, kiut 794 q_s_ b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )794 q_s_1d(ji,jk) = rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) 795 795 END DO 796 796 END DO -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
r4688 r4872 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 -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r4870 r4872 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 76 INTEGER :: ji,jj,jk,jl ! dummy loop indices … … 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(:) :: zv_frazb ! accretion of frazil ice at the ice bottom104 REAL(wp), POINTER, DIMENSION(:) :: zv_frazb ! accretion of frazil ice at the ice bottom 105 105 REAL(wp), POINTER, DIMENSION(:) :: zvrel_1d ! relative ice / frazil velocity (1D vector) 106 106 107 REAL(wp), POINTER, DIMENSION(:,:) :: zv_ old! old volume of ice in category jl108 REAL(wp), POINTER, DIMENSION(:,:) :: za_ old! old area of ice in category jl109 REAL(wp), POINTER, DIMENSION(:,:) :: za_i_1d 110 REAL(wp), POINTER, DIMENSION(:,:) :: zv_i_1d 111 REAL(wp), POINTER, DIMENSION(:,:) :: zoa_i_1d 112 REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_1d 113 114 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 115 115 116 116 REAL(wp), POINTER, DIMENSION(:,:) :: zvrel ! relative ice / frazil velocity … … 120 120 CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 121 121 CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 122 CALL wrk_alloc( jpij,jpl, zv_ old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )122 CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 123 123 CALL wrk_alloc( jpij,jkmax,jpl, ze_i_1d ) 124 124 CALL wrk_alloc( jpi,jpj, zvrel ) … … 302 302 303 303 CALL tab_2d_1d( nbpac, qlead_1d (1:nbpac) , qlead , jpi, jpj, npac(1:nbpac) ) 304 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) ) 305 305 CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac) , sfx_opw, jpi, jpj, npac(1:nbpac) ) 306 306 CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw, jpi, jpj, npac(1:nbpac) ) 307 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) ) 308 308 CALL tab_2d_1d( nbpac, zvrel_1d (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 309 309 … … 318 318 ! Keep old ice areas and volume in memory 319 319 !----------------------------------------- 320 zv_ old(1:nbpac,:) = zv_i_1d(1:nbpac,:)321 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,:) 322 322 !---------------------- 323 323 ! Thickness of new ice … … 326 326 zh_newice(ji) = hiccrit 327 327 END DO 328 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) 329 329 330 330 !---------------------- … … 350 350 DO ji = 1, nbpac 351 351 ztmelts = - tmut * zs_newice(ji) + rtt ! Melting point (K) 352 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_ b(ji) ) &353 & + 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 ) ) & 354 354 & - rcp * ( ztmelts - rtt ) ) 355 355 END DO ! ji … … 369 369 zEi = - ze_newice(ji) / rhoic ! specific enthalpy of forming ice [J/kg] 370 370 371 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] 372 372 ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 373 373 … … 440 440 DO ji = 1, nbpac 441 441 jl = jcat(ji) ! categroy in which new ice is put 442 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 443 443 END DO 444 444 … … 448 448 zinda = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 449 449 ze_i_1d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + & 450 & ( 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) ) & 451 451 & * zinda / MAX( zv_i_1d(ji,jl), epsi20 ) 452 452 END DO … … 489 489 DO ji = 1, nbpac 490 490 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) ) ! 0 if no ice and 1 if yes 491 zoa_i_1d(ji,jl) = za_ old(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb491 zoa_i_1d(ji,jl) = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb 492 492 END DO 493 493 END DO … … 498 498 DO jl = 1, jpl 499 499 DO ji = 1, nbpac 500 zdv = zv_i_1d(ji,jl) - zv_ old(ji,jl)500 zdv = zv_i_1d(ji,jl) - zv_b(ji,jl) 501 501 zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 502 502 END DO … … 541 541 CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 542 542 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 543 CALL wrk_dealloc( jpij,jpl, zv_ old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )543 CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 544 544 CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_1d ) 545 545 CALL wrk_dealloc( jpi,jpj, zvrel ) -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
r4688 r4872 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 -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
r4869 r4872 150 150 ! Diagnostics 151 151 ! ------------------------------------------------- 152 d_u_ice_dyn(:,:) = u_ice(:,:) - old_u_ice(:,:)153 d_v_ice_dyn(:,:) = v_ice(:,:) - old_v_ice(:,:)154 d_a_i_trp (:,:,:) = a_i (:,:,:) - old_a_i(:,:,:)155 d_v_s_trp (:,:,:) = v_s (:,:,:) - old_v_s(:,:,:)156 d_v_i_trp (:,:,:) = v_i (:,:,:) - old_v_i(:,:,:)157 d_e_s_trp (:,:,:,:) = e_s (:,:,:,:) - old_e_s(:,:,:,:)158 d_e_i_trp (:,:,1:nlay_i,:) = e_i (:,:,1:nlay_i,:) - old_e_i(:,:,1:nlay_i,:)159 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 (:,:,:) 160 160 d_smv_i_trp(:,:,:) = 0._wp 161 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(:,:,:) 162 162 163 163 ! conservation test … … 175 175 CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_update1 : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :') 176 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 :') 177 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 :') 178 178 179 179 DO jl = 1, jpl … … 188 188 CALL prt_ctl(tab2d_1=o_i (:,:,jl) , clinfo1= ' lim_update1 : o_i : ') 189 189 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_update1 : a_i : ') 190 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 : ') 191 191 CALL prt_ctl(tab2d_1=d_a_i_trp (:,:,jl) , clinfo1= ' lim_update1 : d_a_i_trp : ') 192 192 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_update1 : v_i : ') 193 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 : ') 194 194 CALL prt_ctl(tab2d_1=d_v_i_trp (:,:,jl) , clinfo1= ' lim_update1 : d_v_i_trp : ') 195 195 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_update1 : v_s : ') 196 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 : ') 197 197 CALL prt_ctl(tab2d_1=d_v_s_trp (:,:,jl) , clinfo1= ' lim_update1 : d_v_s_trp : ') 198 198 CALL prt_ctl(tab2d_1=e_i (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1 : e_i1 : ') 199 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 : ') 200 200 CALL prt_ctl(tab2d_1=d_e_i_trp (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1 : de_i1_trp : ') 201 201 CALL prt_ctl(tab2d_1=e_i (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1 : e_i2 : ') 202 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 : ') 203 203 CALL prt_ctl(tab2d_1=d_e_i_trp (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1 : de_i2_trp : ') 204 204 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_update1 : e_snow : ') 205 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 : ') 206 206 CALL prt_ctl(tab2d_1=d_e_s_trp (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1 : d_e_s_trp : ') 207 207 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_update1 : smv_i : ') 208 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 : ') 209 209 CALL prt_ctl(tab2d_1=d_smv_i_trp(:,:,jl) , clinfo1= ' lim_update1 : d_smv_i_trp : ') 210 210 CALL prt_ctl(tab2d_1=oa_i (:,:,jl) , clinfo1= ' lim_update1 : oa_i : ') 211 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 : ') 212 212 CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl) , clinfo1= ' lim_update1 : d_oa_i_trp : ') 213 213 -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
r4869 r4872 182 182 ! Diagnostics 183 183 ! ------------------------------------------------- 184 d_a_i_thd(:,:,:) = a_i(:,:,:) - old_a_i(:,:,:)185 d_v_s_thd(:,:,:) = v_s(:,:,:) - old_v_s(:,:,:)186 d_v_i_thd(:,:,:) = v_i(:,:,:) - old_v_i(:,:,:)187 d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:)188 d_e_i_thd(:,:,1:nlay_i,:) = e_i(:,:,1:nlay_i,:) - old_e_i(:,:,1:nlay_i,:)189 !?? 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 (:,:,:) 190 190 d_smv_i_thd(:,:,:) = 0._wp 191 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(:,:,:) 192 192 ! diag only (clem) 193 193 dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday … … 215 215 CALL prt_ctl(tab2d_1=strength , clinfo1=' lim_update2 : strength :') 216 216 CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_update2 : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :') 217 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 :') 218 218 219 219 DO jl = 1, jpl … … 228 228 CALL prt_ctl(tab2d_1=o_i (:,:,jl) , clinfo1= ' lim_update2 : o_i : ') 229 229 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_update2 : a_i : ') 230 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 : ') 231 231 CALL prt_ctl(tab2d_1=d_a_i_thd (:,:,jl) , clinfo1= ' lim_update2 : d_a_i_thd : ') 232 232 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_update2 : v_i : ') 233 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 : ') 234 234 CALL prt_ctl(tab2d_1=d_v_i_thd (:,:,jl) , clinfo1= ' lim_update2 : d_v_i_thd : ') 235 235 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_update2 : v_s : ') 236 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 : ') 237 237 CALL prt_ctl(tab2d_1=d_v_s_thd (:,:,jl) , clinfo1= ' lim_update2 : d_v_s_thd : ') 238 238 CALL prt_ctl(tab2d_1=e_i (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2 : e_i1 : ') 239 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 : ') 240 240 CALL prt_ctl(tab2d_1=d_e_i_thd (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2 : de_i1_thd : ') 241 241 CALL prt_ctl(tab2d_1=e_i (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2 : e_i2 : ') 242 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 : ') 243 243 CALL prt_ctl(tab2d_1=d_e_i_thd (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2 : de_i2_thd : ') 244 244 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_update2 : e_snow : ') 245 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 : ') 246 246 CALL prt_ctl(tab2d_1=d_e_s_thd (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2 : d_e_s_thd : ') 247 247 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_update2 : smv_i : ') 248 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 : ') 249 249 CALL prt_ctl(tab2d_1=d_smv_i_thd(:,:,jl) , clinfo1= ' lim_update2 : d_smv_i_thd : ') 250 250 CALL prt_ctl(tab2d_1=oa_i (:,:,jl) , clinfo1= ' lim_update2 : oa_i : ') 251 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 : ') 252 252 CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl) , clinfo1= ' lim_update2 : d_oa_i_thd : ') 253 253 -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r4688 r4872 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 -
trunk/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
r4866 r4872 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 … … 86 86 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sprecip_1d !: <==> the 2D sprecip 87 87 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: frld_1d !: <==> the 2D frld 88 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: at_i_ b!: <==> the 2D at_i88 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: at_i_1d !: <==> the 2D at_i 89 89 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fhtur_1d !: <==> the 2D fhtur 90 90 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fhld_1d !: <==> the 2D fhld … … 99 99 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dsm_i_se_1d !: Ice salinity variations due to basal salt entrapment 100 100 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dsm_i_si_1d !: Ice salinity variations due to lateral accretion 101 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hicol_ b!: Ice collection thickness accumulated in leads101 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hicol_1d !: Ice collection thickness accumulated in leads 102 102 103 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: t_su_ b!: <==> the 2D t_su104 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_i_ b!: <==> the 2D a_i105 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_i_ b!: <==> the 2D ht_s106 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_s_ b!: <==> the 2D ht_i107 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fc_su !: Surface Conduction flux108 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fc_bo_i !: Bottom Conduction flux109 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_ b!: Ice bulk salinity [ppt]114 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 115 115 116 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: iatte_1d !: clemattenuation coef of the input solar flux (unitless)117 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: oatte_1d !: clemattenuation coef of the input solar flux (unitless)116 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: iatte_1d !: attenuation coef of the input solar flux (unitless) 117 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: oatte_1d !: attenuation coef of the input solar flux (unitless) 118 118 119 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_s_ b!: corresponding to the 2D var t_s120 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_i_ b!: corresponding to the 2D var t_i121 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: s_i_ b!: profiled ice salinity122 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_i_ b!: Ice enthalpy per unit volume123 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_s_ b!: Snow enthalpy per unit volume119 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_s_1d !: corresponding to the 2D var t_s 120 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_i_1d !: corresponding to the 2D var t_i 121 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: s_i_1d !: profiled ice salinity 122 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_i_1d !: Ice enthalpy per unit volume 123 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_s_1d !: Snow enthalpy per unit volume 124 124 125 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qh_i_old 126 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_i_old 125 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qh_i_old !: ice heat content (q*h, J.m-2) 126 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_i_old !: ice thickness layer (m) 127 127 128 128 INTEGER , PUBLIC :: jiindex_1d ! 1D index of debugging point … … 148 148 & qsr_ice_1d (jpij) , & 149 149 & fr1_i0_1d(jpij) , fr2_i0_1d(jpij) , qns_ice_1d(jpij) , & 150 & t_bo_ b(jpij) , iatte_1d (jpij) , oatte_1d (jpij) , &150 & t_bo_1d (jpij) , iatte_1d (jpij) , oatte_1d (jpij) , & 151 151 & hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij), & 152 152 & hfx_dif_1d(jpij) ,hfx_opw_1d(jpij) , & … … 155 155 & hfx_res_1d(jpij) , hfx_err_rem_1d(jpij), STAT=ierr(1) ) 156 156 ! 157 ALLOCATE( sprecip_1d (jpij) , frld_1d (jpij) , at_i_ b(jpij) , &157 ALLOCATE( sprecip_1d (jpij) , frld_1d (jpij) , at_i_1d (jpij) , & 158 158 & fhtur_1d (jpij) , wfx_snw_1d (jpij) , wfx_spr_1d (jpij) , & 159 159 & fhld_1d (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d(jpij) , wfx_bom_1d(jpij) , & … … 165 165 & sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , & 166 166 & dsm_i_fl_1d(jpij) , dsm_i_gd_1d(jpij) , dsm_i_se_1d(jpij) , & 167 & dsm_i_si_1d(jpij) , hicol_ b(jpij) , STAT=ierr(2) )167 & dsm_i_si_1d(jpij) , hicol_1d (jpij) , STAT=ierr(2) ) 168 168 ! 169 ALLOCATE( t_su_ b (jpij) , a_i_b (jpij) , ht_i_b(jpij) , &170 & ht_s_ b(jpij) , fc_su (jpij) , fc_bo_i (jpij) , &169 ALLOCATE( t_su_1d (jpij) , a_i_1d (jpij) , ht_i_1d (jpij) , & 170 & ht_s_1d (jpij) , fc_su (jpij) , fc_bo_i (jpij) , & 171 171 & dh_s_tot (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) , & 172 & dh_snowice(jpij) , sm_i_ b(jpij) , s_i_new (jpij) , &173 & t_s_ b(jpij,nlay_s), &174 & t_i_ b(jpij,jkmax), s_i_b(jpij,jkmax) , &175 & q_i_ b(jpij,jkmax), q_s_b(jpij,jkmax) , &172 & dh_snowice(jpij) , sm_i_1d (jpij) , s_i_new (jpij) , & 173 & t_s_1d(jpij,nlay_s), & 174 & t_i_1d(jpij,jkmax), s_i_1d(jpij,jkmax) , & 175 & q_i_1d(jpij,jkmax), q_s_1d(jpij,jkmax) , & 176 176 & qh_i_old(jpij,0:jkmax), h_i_old(jpij,0:jkmax) , STAT=ierr(3)) 177 177 !
Note: See TracChangeset
for help on using the changeset viewer.
