- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icevar.F90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/ICE/icevar.F90
r12178 r12928 82 82 END INTERFACE 83 83 84 !! * Substitutions 85 # include "do_loop_substitute.h90" 84 86 !!---------------------------------------------------------------------- 85 87 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 115 117 ! 116 118 ato_i(:,:) = 1._wp - at_i(:,:) ! open water fraction 117 119 ! 120 !!GS: tm_su always needed by ABL over sea-ice 121 ALLOCATE( z1_at_i(jpi,jpj) ) 122 WHERE( at_i(:,:) > epsi20 ) ; z1_at_i(:,:) = 1._wp / at_i(:,:) 123 ELSEWHERE ; z1_at_i(:,:) = 0._wp 124 END WHERE 125 tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 126 WHERE( at_i(:,:)<=epsi20 ) tm_su(:,:) = rt0 127 ! 118 128 ! The following fields are calculated for diagnostics and outputs only 119 129 ! ==> Do not use them for other purposes 120 130 IF( kn > 1 ) THEN 121 131 ! 122 ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) ) 123 WHERE( at_i(:,:) > epsi20 ) ; z1_at_i(:,:) = 1._wp / at_i(:,:) 124 ELSEWHERE ; z1_at_i(:,:) = 0._wp 125 END WHERE 132 ALLOCATE( z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) ) 126 133 WHERE( vt_i(:,:) > epsi20 ) ; z1_vt_i(:,:) = 1._wp / vt_i(:,:) 127 134 ELSEWHERE ; z1_vt_i(:,:) = 0._wp … … 136 143 ! 137 144 ! ! mean temperature (K), salinity and age 138 tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)139 145 tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 140 146 om_i (:,:) = SUM( oa_i(:,:,:) , dim=3 ) * z1_at_i(:,:) … … 154 160 ! ! put rt0 where there is no ice 155 161 WHERE( at_i(:,:)<=epsi20 ) 156 tm_su(:,:) = rt0157 162 tm_si(:,:) = rt0 158 163 tm_i (:,:) = rt0 … … 165 170 END WHERE 166 171 ! 167 DEALLOCATE( z1_ at_i , z1_vt_i , z1_vt_s )172 DEALLOCATE( z1_vt_i , z1_vt_s ) 168 173 ! 169 174 ENDIF 175 ! 176 DEALLOCATE( z1_at_i ) 170 177 ! 171 178 END SUBROUTINE ice_var_agg … … 236 243 zlay_i = REAL( nlay_i , wp ) ! number of layers 237 244 DO jl = 1, jpl 238 DO jk = 1, nlay_i 239 DO jj = 1, jpj 240 DO ji = 1, jpi 241 IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area 242 ! 243 ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] 244 ztmelts = - sz_i(ji,jj,jk,jl) * rTmlt ! Ice layer melt temperature [C] 245 ! Conversion q(S,T) -> T (second order equation) 246 zbbb = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus 247 zccc = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) ) 248 t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts 249 ! 250 ELSE !--- no ice 251 t_i(ji,jj,jk,jl) = rt0 252 ENDIF 253 END DO 254 END DO 255 END DO 245 DO_3D_11_11( 1, nlay_i ) 246 IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area 247 ! 248 ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] 249 ztmelts = - sz_i(ji,jj,jk,jl) * rTmlt ! Ice layer melt temperature [C] 250 ! Conversion q(S,T) -> T (second order equation) 251 zbbb = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus 252 zccc = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) ) 253 t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts 254 ! 255 ELSE !--- no ice 256 t_i(ji,jj,jk,jl) = rt0 257 ENDIF 258 END_3D 256 259 END DO 257 260 … … 344 347 z1_dS = 1._wp / ( zsi1 - zsi0 ) 345 348 DO jl = 1, jpl 346 DO jj = 1, jpj 347 DO ji = 1, jpi 348 zalpha(ji,jj,jl) = MAX( 0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp ) ) 349 ! ! force a constant profile when SSS too low (Baltic Sea) 350 IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp 351 END DO 352 END DO 349 DO_2D_11_11 350 zalpha(ji,jj,jl) = MAX( 0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp ) ) 351 ! ! force a constant profile when SSS too low (Baltic Sea) 352 IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp 353 END_2D 353 354 END DO 354 355 ! 355 356 ! Computation of the profile 356 357 DO jl = 1, jpl 357 DO jk = 1, nlay_i 358 DO jj = 1, jpj 359 DO ji = 1, jpi 360 ! ! linear profile with 0 surface value 361 zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i 362 zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl) ! weighting the profile 363 sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) ) 364 END DO 365 END DO 366 END DO 358 DO_3D_11_11( 1, nlay_i ) 359 ! ! linear profile with 0 surface value 360 zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i 361 zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl) ! weighting the profile 362 sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) ) 363 END_3D 367 364 END DO 368 365 ! … … 489 486 ! Zap ice energy and use ocean heat to melt ice 490 487 !----------------------------------------------------------------- 491 DO jk = 1, nlay_i 492 DO jj = 1 , jpj 493 DO ji = 1 , jpi 494 ! update exchanges with ocean 495 hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0 496 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj) 497 t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 498 END DO 499 END DO 500 END DO 501 ! 502 DO jk = 1, nlay_s 503 DO jj = 1 , jpj 504 DO ji = 1 , jpi 505 ! update exchanges with ocean 506 hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0 507 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj) 508 t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 509 END DO 510 END DO 511 END DO 488 DO_3D_11_11( 1, nlay_i ) 489 ! update exchanges with ocean 490 hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0 491 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj) 492 t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 493 END_3D 494 ! 495 DO_3D_11_11( 1, nlay_s ) 496 ! update exchanges with ocean 497 hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0 498 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj) 499 t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 500 END_3D 512 501 ! 513 502 !----------------------------------------------------------------- 514 503 ! zap ice and snow volume, add water and salt to ocean 515 504 !----------------------------------------------------------------- 516 DO jj = 1 , jpj 517 DO ji = 1 , jpi 518 ! update exchanges with ocean 519 sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoi * r1_rdtice 520 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoi * r1_rdtice 521 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhos * r1_rdtice 522 ! 523 a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj) 524 v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj) 525 v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj) 526 t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) ) 527 oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj) 528 sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj) 529 ! 530 h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj) 531 h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj) 532 ! 533 a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 534 v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 535 ! 536 END DO 537 END DO 505 DO_2D_11_11 506 ! update exchanges with ocean 507 sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoi * r1_Dt_ice 508 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoi * r1_Dt_ice 509 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhos * r1_Dt_ice 510 ! 511 a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj) 512 v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj) 513 v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj) 514 t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) ) 515 oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj) 516 sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj) 517 ! 518 h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj) 519 h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj) 520 ! 521 a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 522 v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 523 ! 524 END_2D 538 525 ! 539 526 END DO … … 587 574 ! zap ice energy and send it to the ocean 588 575 !---------------------------------------- 589 DO jk = 1, nlay_i 590 DO jj = 1 , jpj 591 DO ji = 1 , jpi 592 IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 593 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0 594 pe_i(ji,jj,jk,jl) = 0._wp 595 ENDIF 596 END DO 597 END DO 598 END DO 599 ! 600 DO jk = 1, nlay_s 601 DO jj = 1 , jpj 602 DO ji = 1 , jpi 603 IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 604 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0 605 pe_s(ji,jj,jk,jl) = 0._wp 606 ENDIF 607 END DO 608 END DO 609 END DO 576 DO_3D_11_11( 1, nlay_i ) 577 IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 578 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0 579 pe_i(ji,jj,jk,jl) = 0._wp 580 ENDIF 581 END_3D 582 ! 583 DO_3D_11_11( 1, nlay_s ) 584 IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 585 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0 586 pe_s(ji,jj,jk,jl) = 0._wp 587 ENDIF 588 END_3D 610 589 ! 611 590 !----------------------------------------------------- 612 591 ! zap ice and snow volume, add water and salt to ocean 613 592 !----------------------------------------------------- 614 DO jj = 1 , jpj 615 DO ji = 1 , jpi 616 IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 617 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt 618 pv_i (ji,jj,jl) = 0._wp 619 ENDIF 620 IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 621 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt 622 pv_s (ji,jj,jl) = 0._wp 623 ENDIF 624 IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 625 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt 626 psv_i (ji,jj,jl) = 0._wp 627 ENDIF 628 END DO 629 END DO 593 DO_2D_11_11 594 IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 595 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt 596 pv_i (ji,jj,jl) = 0._wp 597 ENDIF 598 IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 599 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt 600 pv_s (ji,jj,jl) = 0._wp 601 ENDIF 602 IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp .OR. pv_i(ji,jj,jl) <= 0._wp ) THEN 603 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt 604 psv_i (ji,jj,jl) = 0._wp 605 ENDIF 606 END_2D 630 607 ! 631 608 END DO … … 740 717 !! ** Purpose : compute the equivalent ssh in lead when sea ice is embedded 741 718 !! 742 !! ** Method : ssh_lead = ssh + (Mice + Msnow) / r au0719 !! ** Method : ssh_lead = ssh + (Mice + Msnow) / rho0 743 720 !! 744 721 !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira, … … 770 747 zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 771 748 ! 772 zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_r au0749 zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rho0 773 750 ! 774 751 ELSE … … 960 937 ! In case snow load is in excess that would lead to transformation from snow to ice 961 938 ! Then, transfer the snow excess into the ice (different from icethd_dh) 962 zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - r au0 ) * ph_i(ji,jl) ) * r1_rau0 )939 zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rho0 ) * ph_i(ji,jl) ) * r1_rho0 ) 963 940 ! recompute h_i, h_s avoiding out of bounds values 964 941 ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh )
Note: See TracChangeset
for help on using the changeset viewer.
