Changeset 8325 for branches/2017/dev_r8183_ICEMODEL/NEMOGCM
- Timestamp:
- 2017-07-12T16:51:20+02:00 (7 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 9 edited
-
ice.F90 (modified) (1 diff)
-
limthd.F90 (modified) (7 diffs)
-
limthd_da.F90 (modified) (1 diff)
-
limthd_dh.F90 (modified) (22 diffs)
-
limthd_dif.F90 (modified) (5 diffs)
-
limthd_ent.F90 (modified) (6 diffs)
-
limthd_lac.F90 (modified) (2 diffs)
-
limvar.F90 (modified) (3 diffs)
-
thd_ice.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r8321 r8325 110 110 !! oa_i ! - ! Sea ice areal age content | s | 111 111 !! e_i ! - ! Ice enthalpy | J/m2 | 112 !! - ! q_i_1d ! Ice enthalpy per unit vol. | J/m3 |112 !! - ! e_i_1d ! Ice enthalpy per unit vol. | J/m3 | 113 113 !! e_s ! - ! Snow enthalpy | J/m2 | 114 !! - ! q_s_1d ! Snow enthalpy per unit vol. | J/m3 |114 !! - ! e_s_1d ! Snow enthalpy per unit vol. | J/m3 | 115 115 !! | 116 116 !!-------------|-------------|---------------------------------|-------| -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r8321 r8325 131 131 CALL lbc_lnk( zfric, 'T', 1. ) 132 132 ! 133 !----------------------------------! 134 ! Initialization and units change 135 !----------------------------------! 136 ftr_ice(:,:,:) = 0._wp ! part of solar radiation transmitted through the ice 137 138 ! Change the units of heat content; from J/m2 to J/m3 139 DO jl = 1, jpl 140 DO jk = 1, nlay_i 141 DO jj = 1, jpj 142 DO ji = 1, jpi 143 rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 ) ) 144 !Energy of melting q(S,T) [J.m-3] 145 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i ) 146 END DO 147 END DO 148 END DO 149 DO jk = 1, nlay_s 150 DO jj = 1, jpj 151 DO ji = 1, jpi 152 rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) ) 153 !Energy of melting q(S,T) [J.m-3] 154 e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s ) 155 END DO 156 END DO 157 END DO 158 END DO 133 ftr_ice(:,:,:) = 0._wp ! initialization (part of solar radiation transmitted through the ice) 159 134 160 135 !--------------------------------------------------------------------! … … 269 244 CALL lim_thd_1d2d( nbpb, jl, 1 ) ! --- Move to 1D arrays --- ! 270 245 ! 246 DO jk = 1, nlay_i ! --- Change units from J/m2 to J/m3 --- ! 247 WHERE( ht_i_1d(:)>0._wp ) e_i_1d(:,jk) = e_i_1d(:,jk) / (ht_i_1d(:) * a_i_1d(:)) * nlay_i 248 ENDDO 249 DO jk = 1, nlay_s 250 WHERE( ht_s_1d(:)>0._wp ) e_s_1d(:,jk) = e_s_1d(:,jk) / (ht_s_1d(:) * a_i_1d(:)) * nlay_s 251 ENDDO 252 ! 271 253 IF( ln_limdH ) CALL lim_thd_dif( 1, nbpb ) ! --- Ice/Snow Temperature profile --- ! 272 254 ! 273 255 IF( ln_limdH ) CALL lim_thd_dh( 1, nbpb ) ! --- Ice/Snow thickness --- ! 274 256 ! 275 IF( ln_limdH ) CALL lim_thd_ent( 1, nbpb, q_i_1d(1:nbpb,:) ) ! --- Ice enthalpy remapping --- !257 IF( ln_limdH ) CALL lim_thd_ent( 1, nbpb, e_i_1d(1:nbpb,:) ) ! --- Ice enthalpy remapping --- ! 276 258 ! 277 259 CALL lim_thd_sal( 1, nbpb ) ! --- Ice salinity --- ! … … 285 267 END IF 286 268 ! 269 DO jk = 1, nlay_i ! --- Change units from J/m3 to J/m2 --- ! 270 e_i_1d(:,jk) = e_i_1d(:,jk) * ht_i_1d(:) * a_i_1d(:) * r1_nlay_i 271 ENDDO 272 DO jk = 1, nlay_s 273 e_s_1d(:,jk) = e_s_1d(:,jk) * ht_s_1d(:) * a_i_1d(:) * r1_nlay_s 274 ENDDO 275 ! 287 276 CALL lim_thd_1d2d( nbpb, jl, 2 ) ! --- Move to 2D arrays --- ! 288 277 ! … … 294 283 IF( ln_limdA) CALL lim_thd_da ! --- lateral melting --- ! 295 284 296 ! Enthalpies are global variables we have to readjust the units (heat content in J/m2)297 DO jl = 1, jpl298 DO jk = 1, nlay_i299 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) * r1_nlay_i300 END DO301 DO jk = 1, nlay_s302 e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) * r1_nlay_s303 END DO304 END DO305 306 285 ! Change thickness to volume 307 286 v_i(:,:,:) = ht_i(:,:,:) * a_i(:,:,:) … … 375 354 ! Conversion q(S,T) -> T (second order equation) 376 355 zaaa = cpic 377 zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + q_i_1d(ji,jk) * r1_rhoic - lfus356 zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + e_i_1d(ji,jk) * r1_rhoic - lfus 378 357 zccc = lfus * ( ztmelts - rt0 ) 379 358 zdiscrim = SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) … … 451 430 DO jk = 1, nlay_s 452 431 CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 453 CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )432 CALL tab_2d_1d( nbpb, e_s_1d(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 454 433 END DO 455 434 DO jk = 1, nlay_i 456 435 CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 457 CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )436 CALL tab_2d_1d( nbpb, e_i_1d(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 458 437 CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) ) 459 438 END DO … … 523 502 DO jk = 1, nlay_s 524 503 CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d (1:nbpb,jk), jpi, jpj) 525 CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d (1:nbpb,jk), jpi, jpj)504 CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, e_s_1d (1:nbpb,jk), jpi, jpj) 526 505 END DO 527 506 DO jk = 1, nlay_i 528 507 CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d (1:nbpb,jk), jpi, jpj) 529 CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d (1:nbpb,jk), jpi, jpj)508 CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, e_i_1d (1:nbpb,jk), jpi, jpj) 530 509 CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d (1:nbpb,jk), jpi, jpj) 531 510 END DO -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_da.F90
r7753 r8325 150 150 151 151 ! Contribution to heat flux into the ocean [W.m-2], <0 152 hfx_thd(ji,jj) = hfx_thd(ji,jj) + zda * r1_rdtice * ( ht_i(ji,jj,jl) * SUM( e_i(ji,jj,:,jl) ) * r1_nlay_i & 153 & + ht_s(ji,jj,jl) * e_s(ji,jj,1,jl) * r1_nlay_s ) 152 !clemX hfx_thd(ji,jj) = hfx_thd(ji,jj) + zda * r1_rdtice * ( ht_i(ji,jj,jl) * SUM( e_i(ji,jj,:,jl) ) * r1_nlay_i & 153 ! & + ht_s(ji,jj,jl) * e_s(ji,jj,1,jl) * r1_nlay_s ) 154 hfx_thd(ji,jj) = hfx_thd(ji,jj) + rswitch * zda_tot(ji,jj) / MAX( at_i(ji,jj), epsi10 ) & 155 & * r1_rdtice * ( SUM( e_i(ji,jj,:,jl) ) + e_s(ji,jj,1,jl) ) 154 156 155 157 ! Contribution to mass flux -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r8316 r8325 104 104 INTEGER , POINTER, DIMENSION(:,:) :: icount ! number of layers vanished by melting 105 105 106 REAL(wp), POINTER, DIMENSION(:) :: z qh_i ! total ice heat content (J.m-2)106 REAL(wp), POINTER, DIMENSION(:) :: zeh_i ! total ice heat content (J.m-2) 107 107 REAL(wp), POINTER, DIMENSION(:) :: zsnw ! distribution of snow after wind blowing 108 108 … … 121 121 122 122 CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema ) 123 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, z qh_i )123 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zeh_i ) 124 124 CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 125 125 CALL wrk_alloc( jpij, nlay_i, icount ) … … 127 127 zqprec (:) = 0._wp ; zq_su (:) = 0._wp ; zq_bo (:) = 0._wp ; zf_tt(:) = 0._wp 128 128 zq_rema (:) = 0._wp ; zsnw (:) = 0._wp ; zevap_rema(:) = 0._wp ; 129 zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; z qh_i(:) = 0._wp129 zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; zeh_i(:) = 0._wp 130 130 131 131 zdeltah(:,:) = 0._wp ; zh_i(:,:) = 0._wp … … 134 134 ! Initialize enthalpy at nlay_i+1 135 135 DO ji = kideb, kiut 136 q_i_1d(ji,nlay_i+1) = 0._wp136 e_i_1d(ji,nlay_i+1) = 0._wp 137 137 END DO 138 138 139 139 ! initialize layer thicknesses and enthalpies 140 140 h_i_old (:,0:nlay_i+1) = 0._wp 141 qh_i_old(:,0:nlay_i+1) = 0._wp141 eh_i_old(:,0:nlay_i+1) = 0._wp 142 142 DO jk = 1, nlay_i 143 143 DO ji = kideb, kiut 144 144 h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i 145 qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk)145 eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk) 146 146 ENDDO 147 147 ENDDO … … 167 167 IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting 168 168 ! Contribution to heat flux to the ocean [W.m-2], < 0 169 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice169 hfx_res_1d(ji) = hfx_res_1d(ji) + e_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 170 170 ! Contribution to mass flux 171 171 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 172 172 ! updates 173 173 ht_s_1d(ji) = 0._wp 174 q_s_1d (ji,1) = 0._wp174 e_s_1d (ji,1) = 0._wp 175 175 t_s_1d (ji,1) = rt0 176 176 END IF … … 184 184 DO ji = kideb, kiut 185 185 zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i 186 z qh_i(ji) = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk)186 zeh_i(ji) = zeh_i(ji) + e_i_1d(ji,jk) * zh_i(ji,jk) 187 187 END DO 188 188 END DO … … 248 248 ! thickness change 249 249 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 250 rswitch = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) )251 zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 )250 rswitch = rswitch * ( MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) ) 251 zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( e_s_1d(ji,jk), epsi20 ) 252 252 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting 253 253 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 254 254 ! heat used to melt snow(W.m-2, >0) 255 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice255 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d(ji,jk) * r1_rdtice 256 256 ! snow melting only = water into the ocean (then without snow precip) 257 257 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 258 258 ! updates available heat + thickness 259 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )259 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 260 260 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 261 261 END DO … … 274 274 ! Heat flux by sublimation [W.m-2], < 0 (sublimate first snow that had fallen, then pre-existing snow) 275 275 zdeltah(ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 276 hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1) &276 hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) & 277 277 & ) * a_i_1d(ji) * r1_rdtice 278 278 ! Mass flux by sublimation … … 298 298 DO ji = kideb,kiut 299 299 rswitch = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) ) 300 q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) * &300 e_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) * & 301 301 & ( ( zdh_s_pre(ji) ) * zqprec(ji) + & 302 302 & ( ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) … … 314 314 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 315 315 316 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0]316 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 317 317 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 318 318 ! set up at 0 since no energy is needed to melt water...(it is already melted) … … 336 336 ELSE !!! Surface melting 337 337 338 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0]338 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 339 339 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] 340 340 zdE = zEi - zEw ! Specific enthalpy difference < 0 … … 377 377 sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * sm_i_1d(ji) * r1_rdtice 378 378 ! Heat flux [W.m-2], < 0 379 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * q_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice379 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice 380 380 ! Mass flux > 0 381 381 wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice … … 392 392 393 393 ! update heat content (J.m-2) and layer thickness 394 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)394 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 395 395 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 396 396 END DO … … 464 464 dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 465 465 466 q_i_1d(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0)466 e_i_1d(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0) 467 467 468 468 ENDIF … … 502 502 503 503 ! update heat content (J.m-2) and layer thickness 504 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)504 eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * e_i_1d(ji,nlay_i+1) 505 505 h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji) 506 506 ENDIF … … 519 519 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 520 520 521 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)521 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 522 522 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 523 523 ! set up at 0 since no energy is needed to melt water...(it is already melted) … … 539 539 540 540 ! update heat content (J.m-2) and layer thickness 541 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)541 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 542 542 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 543 543 544 544 ELSE !!! Basal melting 545 545 546 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)546 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 547 547 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of meltwater (J/kg, <0) 548 548 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) … … 575 575 576 576 ! update heat content (J.m-2) and layer thickness 577 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)577 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 578 578 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 579 579 ENDIF … … 599 599 zq_rema(ji) = zq_su(ji) + zq_bo(ji) 600 600 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) ! =1 if snow 601 rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) )602 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 )601 rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) ) 602 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 ) 603 603 zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 604 604 dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1) 605 605 ht_s_1d (ji) = ht_s_1d(ji) + zdeltah(ji,1) 606 606 607 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * q_s_1d(ji,1) ! update available heat (J.m-2)607 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1) ! update available heat (J.m-2) 608 608 ! heat used to melt snow 609 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * q_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0)609 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0) 610 610 ! Contribution to mass flux 611 611 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice … … 661 661 662 662 ! update heat content (J.m-2) and layer thickness 663 qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw663 eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw 664 664 h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 665 665 … … 679 679 ! mask enthalpy 680 680 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) ) ) 681 q_s_1d(ji,jk) = rswitch * q_s_1d(ji,jk)682 ! recalculate t_s_1d from q_s_1d683 t_s_1d(ji,jk) = rt0 + rswitch * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )681 e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk) 682 ! recalculate t_s_1d from e_s_1d 683 t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 684 684 END DO 685 685 END DO … … 689 689 690 690 CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema ) 691 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, z qh_i )691 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zeh_i ) 692 692 CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 693 693 CALL wrk_dealloc( jpij, nlay_i, icount ) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r8313 r8325 109 109 REAL(wp) :: zraext_s = 10._wp ! extinction coefficient of radiation in the snow 110 110 REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity 111 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered a s 0°C111 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered at 0C 112 112 REAL(wp) :: ztmelt_i ! ice melting temperature 113 113 REAL(wp) :: zhsu … … 181 181 zdq(:) = 0._wp ; zq_ini(:) = 0._wp 182 182 DO ji = kideb, kiut 183 zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + &184 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )183 zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & 184 & SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 185 185 END DO 186 186 … … 782 782 ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 783 783 DO ji = kideb, kiut 784 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + &785 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )784 zdq(ji) = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & 785 & SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 786 786 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 787 787 zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice … … 840 840 t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point 841 841 ! (sometimes dif scheme produces abnormally high temperatures) 842 q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) &842 e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) & 843 843 & + lfus * ( 1.0 - ( ztmelts-rt0 ) / ( t_i_1d(ji,jk) - rt0 ) ) & 844 844 & - rcp * ( ztmelts-rt0 ) ) … … 847 847 DO jk = 1, nlay_s ! Snow energy of melting 848 848 DO ji = kideb, kiut 849 q_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )849 e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) 850 850 END DO 851 851 END DO -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
r8316 r8325 75 75 INTEGER :: jk0, jk1 ! old/new layer indices 76 76 ! 77 REAL(wp), POINTER, DIMENSION(:,:) :: z qh_cum0, zh_cum0 ! old cumulative enthlapies and layers interfaces78 REAL(wp), POINTER, DIMENSION(:,:) :: z qh_cum1, zh_cum1 ! new cumulative enthlapies and layers interfaces77 REAL(wp), POINTER, DIMENSION(:,:) :: zeh_cum0, zh_cum0 ! old cumulative enthlapies and layers interfaces 78 REAL(wp), POINTER, DIMENSION(:,:) :: zeh_cum1, zh_cum1 ! new cumulative enthlapies and layers interfaces 79 79 REAL(wp), POINTER, DIMENSION(:) :: zhnew ! new layers thicknesses 80 80 !!------------------------------------------------------------------- 81 81 82 CALL wrk_alloc( jpij, nlay_i+3, z qh_cum0, zh_cum0, kjstart = 0 )83 CALL wrk_alloc( jpij, nlay_i+1, z qh_cum1, zh_cum1, kjstart = 0 )82 CALL wrk_alloc( jpij, nlay_i+3, zeh_cum0, zh_cum0, kjstart = 0 ) 83 CALL wrk_alloc( jpij, nlay_i+1, zeh_cum1, zh_cum1, kjstart = 0 ) 84 84 CALL wrk_alloc( jpij, zhnew ) 85 85 … … 87 87 ! 1) Cumulative integral of old enthalpy * thickness and layers interfaces 88 88 !-------------------------------------------------------------------------- 89 z qh_cum0(:,0:nlay_i+2) = 0._wp89 zeh_cum0(:,0:nlay_i+2) = 0._wp 90 90 zh_cum0 (:,0:nlay_i+2) = 0._wp 91 91 DO jk0 = 1, nlay_i+2 92 92 DO ji = kideb, kiut 93 z qh_cum0(ji,jk0) = zqh_cum0(ji,jk0-1) + qh_i_old(ji,jk0-1)93 zeh_cum0(ji,jk0) = zeh_cum0(ji,jk0-1) + eh_i_old(ji,jk0-1) 94 94 zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + h_i_old (ji,jk0-1) 95 95 ENDDO … … 112 112 ENDDO 113 113 114 z qh_cum1(:,0:nlay_i) = 0._wp114 zeh_cum1(:,0:nlay_i) = 0._wp 115 115 ! new cumulative q*h => linear interpolation 116 116 DO jk0 = 1, nlay_i+1 … … 118 118 DO ji = kideb, kiut 119 119 IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN 120 z qh_cum1(ji,jk1) = ( zqh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1 ) ) + &121 & z qh_cum0(ji,jk0 ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) ) &120 zeh_cum1(ji,jk1) = ( zeh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1 ) ) + & 121 & zeh_cum0(ji,jk0 ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) ) & 122 122 & / ( zh_cum0(ji,jk0) - zh_cum0(ji,jk0-1) ) 123 123 ENDIF … … 126 126 ENDDO 127 127 ! to ensure that total heat content is strictly conserved, set: 128 z qh_cum1(:,nlay_i) = zqh_cum0(:,nlay_i+2)128 zeh_cum1(:,nlay_i) = zeh_cum0(:,nlay_i+2) 129 129 130 130 ! new enthalpies … … 132 132 DO ji = kideb, kiut 133 133 rswitch = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) ) 134 qnew(ji,jk1) = rswitch * ( z qh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 )134 qnew(ji,jk1) = rswitch * ( zeh_cum1(ji,jk1) - zeh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 ) 135 135 ENDDO 136 136 ENDDO 137 137 138 138 ! --- diag error on heat remapping --- ! 139 ! comment: if input h_i_old and qh_i_old are already multiplied by a_i (as in limthd_lac),139 ! comment: if input h_i_old and eh_i_old are already multiplied by a_i (as in limthd_lac), 140 140 ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0 141 141 DO ji = kideb, kiut 142 142 hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_rdtice * & 143 & ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( qh_i_old(ji,0:nlay_i+1) ) )143 & ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_i_old(ji,0:nlay_i+1) ) ) 144 144 END DO 145 145 146 146 ! 147 CALL wrk_dealloc( jpij, nlay_i+3, z qh_cum0, zh_cum0, kjstart = 0 )148 CALL wrk_dealloc( jpij, nlay_i+1, z qh_cum1, zh_cum1, kjstart = 0 )147 CALL wrk_dealloc( jpij, nlay_i+3, zeh_cum0, zh_cum0, kjstart = 0 ) 148 CALL wrk_dealloc( jpij, nlay_i+1, zeh_cum1, zh_cum1, kjstart = 0 ) 149 149 CALL wrk_dealloc( jpij, zhnew ) 150 150 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r8316 r8325 446 446 ! for remapping 447 447 h_i_old (1:nbpac,0:nlay_i+1) = 0._wp 448 qh_i_old(1:nbpac,0:nlay_i+1) = 0._wp448 eh_i_old(1:nbpac,0:nlay_i+1) = 0._wp 449 449 DO jk = 1, nlay_i 450 450 DO ji = 1, nbpac 451 451 h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i 452 qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk)452 eh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk) 453 453 END DO 454 454 END DO … … 462 462 ! for remapping 463 463 h_i_old (ji,nlay_i+1) = zv_newfra 464 qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra464 eh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra 465 465 ENDDO 466 466 ! --- Ice enthalpy remapping --- ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r8316 r8325 162 162 !!------------------------------------------------------------------ 163 163 INTEGER :: ji, jj, jk, jl ! dummy loop indices 164 REAL(wp) :: z q_i, zaaa, zbbb, zccc, zdiscrim ! local scalars165 REAL(wp) :: ztmelts, z q_s, zfac1, zfac2 ! - -164 REAL(wp) :: ze_i, zaaa, zbbb, zccc, zdiscrim ! local scalars 165 REAL(wp) :: ztmelts, ze_s, zfac1, zfac2 ! - - 166 166 !!------------------------------------------------------------------ 167 167 … … 218 218 DO jj = 1, jpj 219 219 DO ji = 1, jpi 220 ! ! Energy of melting q(S,T) [J.m-3]220 ! ! Energy of melting e(S,T) [J.m-3] 221 221 rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes 222 z q_i = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp)222 ze_i = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp) 223 223 ztmelts = -tmut * s_i(ji,jj,jk,jl) + rt0 ! Ice layer melt temperature 224 224 ! 225 225 zaaa = cpic ! Conversion q(S,T) -> T (second order equation) 226 zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + z q_i * r1_rhoic - lfus226 zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + ze_i * r1_rhoic - lfus 227 227 zccc = lfus * (ztmelts-rt0) 228 228 zdiscrim = SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) … … 245 245 !Energy of melting q(S,T) [J.m-3] 246 246 rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes 247 z q_s = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp)247 ze_s = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp) 248 248 ! 249 t_s(ji,jj,jk,jl) = rt0 + rswitch * ( - zfac1 * z q_s + zfac2 )249 t_s(ji,jj,jk,jl) = rt0 + rswitch * ( - zfac1 * ze_s + zfac2 ) 250 250 t_s(ji,jj,jk,jl) = MIN( rt0, MAX( rt0 - 100._wp , t_s(ji,jj,jk,jl) ) ) ! -100 < t_s < rt0 251 251 END DO -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
r8313 r8325 111 111 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_i_1d !: corresponding to the 2D var t_i 112 112 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: s_i_1d !: profiled ice salinity 113 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_i_1d !: Ice enthalpy per unit volume114 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_s_1d !: Snow enthalpy per unit volume113 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e_i_1d !: Ice enthalpy per unit volume 114 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e_s_1d !: Snow enthalpy per unit volume 115 115 116 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qh_i_old !: ice heat content (q*h, J.m-2)116 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: eh_i_old !: ice heat content (q*h, J.m-2) 117 117 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_i_old !: ice thickness layer (m) 118 118 … … 172 172 ii = ii + 1 173 173 ALLOCATE( t_s_1d (jpij,nlay_s) , t_i_1d (jpij,nlay_i) , s_i_1d(jpij,nlay_i) , & 174 & q_i_1d (jpij,nlay_i+1) , q_s_1d (jpij,nlay_s) , &175 & qh_i_old(jpij,0:nlay_i+1) , h_i_old(jpij,0:nlay_i+1) , STAT=ierr(ii) )174 & e_i_1d (jpij,nlay_i+1) , e_s_1d (jpij,nlay_s) , & 175 & eh_i_old(jpij,0:nlay_i+1) , h_i_old(jpij,0:nlay_i+1) , STAT=ierr(ii) ) 176 176 ! 177 177 ii = ii + 1
Note: See TracChangeset
for help on using the changeset viewer.
