- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icevar.F90 (modified) (24 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/ICE/icevar.F90
r11229 r13463 32 32 !! - vt_s(jpi,jpj) 33 33 !! - at_i(jpi,jpj) 34 !! - st_i(jpi,jpj) 34 35 !! - et_s(jpi,jpj) total snow heat content 35 36 !! - et_i(jpi,jpj) total ice thermal content … … 46 47 !! ice_var_zapneg : remove negative ice fields 47 48 !! ice_var_roundoff : remove negative values arising from roundoff erros 48 !! ice_var_itd : convert 1-cat to jpl-cat49 !! ice_var_itd2 : convert N-cat to jpl-cat50 49 !! ice_var_bv : brine volume 51 50 !! ice_var_enthalpy : compute ice and snow enthalpies from temperature 52 51 !! ice_var_sshdyn : compute equivalent ssh in lead 52 !! ice_var_itd : convert N-cat to M-cat 53 53 !!---------------------------------------------------------------------- 54 54 USE dom_oce ! ocean space and time domain … … 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) … … 104 106 ! 105 107 ! ! integrated values 106 vt_i(:,:) = SUM( v_i(:,:,:) , dim=3 ) 107 vt_s(:,:) = SUM( v_s(:,:,:) , dim=3 ) 108 at_i(:,:) = SUM( a_i(:,:,:) , dim=3 ) 109 et_s(:,:) = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 ) 110 et_i(:,:) = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 ) 108 vt_i(:,:) = SUM( v_i (:,:,:) , dim=3 ) 109 vt_s(:,:) = SUM( v_s (:,:,:) , dim=3 ) 110 st_i(:,:) = SUM( sv_i(:,:,:) , dim=3 ) 111 at_i(:,:) = SUM( a_i (:,:,:) , dim=3 ) 112 et_s(:,:) = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 ) 113 et_i(:,:) = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 ) 111 114 ! 112 115 at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds … … 114 117 ! 115 118 ato_i(:,:) = 1._wp - at_i(:,:) ! open water fraction 116 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 ! 117 128 ! The following fields are calculated for diagnostics and outputs only 118 129 ! ==> Do not use them for other purposes 119 130 IF( kn > 1 ) THEN 120 131 ! 121 ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) ) 122 WHERE( at_i(:,:) > epsi20 ) ; z1_at_i(:,:) = 1._wp / at_i(:,:) 123 ELSEWHERE ; z1_at_i(:,:) = 0._wp 124 END WHERE 132 ALLOCATE( z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) ) 125 133 WHERE( vt_i(:,:) > epsi20 ) ; z1_vt_i(:,:) = 1._wp / vt_i(:,:) 126 134 ELSEWHERE ; z1_vt_i(:,:) = 0._wp … … 135 143 ! 136 144 ! ! mean temperature (K), salinity and age 137 tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)138 145 tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 139 146 om_i (:,:) = SUM( oa_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 140 sm_i (:,:) = SUM( sv_i(:,:,:) , dim=3 )* z1_vt_i(:,:)147 sm_i (:,:) = st_i(:,:) * z1_vt_i(:,:) 141 148 ! 142 149 tm_i(:,:) = 0._wp … … 153 160 ! ! put rt0 where there is no ice 154 161 WHERE( at_i(:,:)<=epsi20 ) 155 tm_su(:,:) = rt0156 162 tm_si(:,:) = rt0 157 163 tm_i (:,:) = rt0 158 164 tm_s (:,:) = rt0 159 165 END WHERE 160 161 DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s ) 166 ! 167 ! ! mean melt pond depth 168 WHERE( at_ip(:,:) > epsi20 ) ; hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:) 169 ELSEWHERE ; hm_ip(:,:) = 0._wp 170 END WHERE 171 ! 172 DEALLOCATE( z1_vt_i , z1_vt_s ) 173 ! 162 174 ENDIF 175 ! 176 DEALLOCATE( z1_at_i ) 163 177 ! 164 178 END SUBROUTINE ice_var_agg … … 229 243 zlay_i = REAL( nlay_i , wp ) ! number of layers 230 244 DO jl = 1, jpl 231 DO jk = 1, nlay_i 232 DO jj = 1, jpj 233 DO ji = 1, jpi 234 IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area 235 ! 236 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] 237 ztmelts = - sz_i(ji,jj,jk,jl) * rTmlt ! Ice layer melt temperature [C] 238 ! Conversion q(S,T) -> T (second order equation) 239 zbbb = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus 240 zccc = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) ) 241 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 242 ! 243 ELSE !--- no ice 244 t_i(ji,jj,jk,jl) = rt0 245 ENDIF 246 END DO 247 END DO 248 END DO 245 DO_3D( 1, 1, 1, 1, 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 249 259 END DO 250 260 … … 263 273 ! 264 274 ! integrated values 265 vt_i (:,:) = SUM( v_i , dim=3 )266 vt_s (:,:) = SUM( v_s , dim=3 )267 at_i (:,:) = SUM( a_i , dim=3 )275 vt_i (:,:) = SUM( v_i , dim=3 ) 276 vt_s (:,:) = SUM( v_s , dim=3 ) 277 at_i (:,:) = SUM( a_i , dim=3 ) 268 278 ! 269 279 END SUBROUTINE ice_var_glo2eqv … … 337 347 z1_dS = 1._wp / ( zsi1 - zsi0 ) 338 348 DO jl = 1, jpl 339 DO jj = 1, jpj 340 DO ji = 1, jpi 341 zalpha(ji,jj,jl) = MAX( 0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp ) ) 342 ! ! force a constant profile when SSS too low (Baltic Sea) 343 IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp 344 END DO 345 END DO 349 DO_2D( 1, 1, 1, 1 ) 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 346 354 END DO 347 355 ! 348 356 ! Computation of the profile 349 357 DO jl = 1, jpl 350 DO jk = 1, nlay_i 351 DO jj = 1, jpj 352 DO ji = 1, jpi 353 ! ! linear profile with 0 surface value 354 zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i 355 zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl) ! weighting the profile 356 sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) ) 357 END DO 358 END DO 359 END DO 358 DO_3D( 1, 1, 1, 1, 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 360 364 END DO 361 365 ! … … 482 486 ! Zap ice energy and use ocean heat to melt ice 483 487 !----------------------------------------------------------------- 484 DO jk = 1, nlay_i 485 DO jj = 1 , jpj 486 DO ji = 1 , jpi 487 ! update exchanges with ocean 488 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 489 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj) 490 t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 491 END DO 492 END DO 493 END DO 494 ! 495 DO jk = 1, nlay_s 496 DO jj = 1 , jpj 497 DO ji = 1 , jpi 498 ! update exchanges with ocean 499 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 500 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj) 501 t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 502 END DO 503 END DO 504 END DO 488 DO_3D( 1, 1, 1, 1, 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( 1, 1, 1, 1, 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 505 501 ! 506 502 !----------------------------------------------------------------- 507 503 ! zap ice and snow volume, add water and salt to ocean 508 504 !----------------------------------------------------------------- 509 DO jj = 1 , jpj 510 DO ji = 1 , jpi 511 ! update exchanges with ocean 512 sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoi * r1_rdtice 513 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoi * r1_rdtice 514 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhos * r1_rdtice 515 ! 516 a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj) 517 v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj) 518 v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj) 519 t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) ) 520 oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj) 521 sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj) 522 ! 523 h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj) 524 h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj) 525 ! 526 a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 527 v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 528 ! 529 END DO 530 END DO 505 DO_2D( 1, 1, 1, 1 ) 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 531 525 ! 532 526 END DO 533 527 534 528 ! to be sure that at_i is the sum of a_i(jl) 535 at_i (:,:) = SUM( a_i(:,:,:), dim=3 ) 536 vt_i (:,:) = SUM( v_i(:,:,:), dim=3 ) 529 at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) 530 vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) 531 !!clem add? 532 ! vt_s (:,:) = SUM( v_s (:,:,:), dim=3 ) 533 ! st_i (:,:) = SUM( sv_i(:,:,:), dim=3 ) 534 ! et_s(:,:) = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 ) 535 ! et_i(:,:) = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 ) 536 !!clem 537 537 538 538 ! open water = 1 if at_i=0 … … 574 574 ! zap ice energy and send it to the ocean 575 575 !---------------------------------------- 576 DO jk = 1, nlay_i 577 DO jj = 1 , jpj 578 DO ji = 1 , jpi 579 IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 580 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0 581 pe_i(ji,jj,jk,jl) = 0._wp 582 ENDIF 583 END DO 584 END DO 585 END DO 586 ! 587 DO jk = 1, nlay_s 588 DO jj = 1 , jpj 589 DO ji = 1 , jpi 590 IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 591 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0 592 pe_s(ji,jj,jk,jl) = 0._wp 593 ENDIF 594 END DO 595 END DO 596 END DO 576 DO_3D( 1, 1, 1, 1, 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( 1, 1, 1, 1, 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 597 589 ! 598 590 !----------------------------------------------------- 599 591 ! zap ice and snow volume, add water and salt to ocean 600 592 !----------------------------------------------------- 601 DO jj = 1 , jpj 602 DO ji = 1 , jpi 603 IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 604 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt 605 pv_i (ji,jj,jl) = 0._wp 606 ENDIF 607 IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 608 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt 609 pv_s (ji,jj,jl) = 0._wp 610 ENDIF 611 IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 612 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt 613 psv_i (ji,jj,jl) = 0._wp 614 ENDIF 615 END DO 616 END DO 593 DO_2D( 1, 1, 1, 1 ) 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 617 607 ! 618 608 END DO … … 645 635 !!------------------------------------------------------------------- 646 636 ! 647 WHERE( pa_i (1:npti,:) < 0._wp .AND. pa_i (1:npti,:) > -epsi10 ) pa_i (1:npti,:) = 0._wp ! a_i must be >= 0 648 WHERE( pv_i (1:npti,:) < 0._wp .AND. pv_i (1:npti,:) > -epsi10 ) pv_i (1:npti,:) = 0._wp ! v_i must be >= 0 649 WHERE( pv_s (1:npti,:) < 0._wp .AND. pv_s (1:npti,:) > -epsi10 ) pv_s (1:npti,:) = 0._wp ! v_s must be >= 0 650 WHERE( psv_i(1:npti,:) < 0._wp .AND. psv_i(1:npti,:) > -epsi10 ) psv_i(1:npti,:) = 0._wp ! sv_i must be >= 0 651 WHERE( poa_i(1:npti,:) < 0._wp .AND. poa_i(1:npti,:) > -epsi10 ) poa_i(1:npti,:) = 0._wp ! oa_i must be >= 0 652 WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 ) pe_i (1:npti,:,:) = 0._wp ! e_i must be >= 0 653 WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 ) pe_s (1:npti,:,:) = 0._wp ! e_s must be >= 0 654 IF ( ln_pnd_H12 ) THEN 655 WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0 656 WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0 637 638 WHERE( pa_i (1:npti,:) < 0._wp ) pa_i (1:npti,:) = 0._wp ! a_i must be >= 0 639 WHERE( pv_i (1:npti,:) < 0._wp ) pv_i (1:npti,:) = 0._wp ! v_i must be >= 0 640 WHERE( pv_s (1:npti,:) < 0._wp ) pv_s (1:npti,:) = 0._wp ! v_s must be >= 0 641 WHERE( psv_i(1:npti,:) < 0._wp ) psv_i(1:npti,:) = 0._wp ! sv_i must be >= 0 642 WHERE( poa_i(1:npti,:) < 0._wp ) poa_i(1:npti,:) = 0._wp ! oa_i must be >= 0 643 WHERE( pe_i (1:npti,:,:) < 0._wp ) pe_i (1:npti,:,:) = 0._wp ! e_i must be >= 0 644 WHERE( pe_s (1:npti,:,:) < 0._wp ) pe_s (1:npti,:,:) = 0._wp ! e_s must be >= 0 645 IF( ln_pnd_H12 ) THEN 646 WHERE( pa_ip(1:npti,:) < 0._wp ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0 647 WHERE( pv_ip(1:npti,:) < 0._wp ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0 657 648 ENDIF 658 649 ! … … 727 718 !! ** Purpose : compute the equivalent ssh in lead when sea ice is embedded 728 719 !! 729 !! ** Method : ssh_lead = ssh + (Mice + Msnow) / r au0720 !! ** Method : ssh_lead = ssh + (Mice + Msnow) / rho0 730 721 !! 731 722 !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira, … … 757 748 zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 758 749 ! 759 zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_r au0750 zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rho0 760 751 ! 761 752 ELSE … … 773 764 !! ** Purpose : converting N-cat ice to jpl ice categories 774 765 !!------------------------------------------------------------------- 775 SUBROUTINE ice_var_itd_1c1c( zhti, zhts, zati, zh_i, zh_s, za_i ) 766 SUBROUTINE ice_var_itd_1c1c( phti, phts, pati , ph_i, ph_s, pa_i, & 767 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 776 768 !!------------------------------------------------------------------- 777 769 !! ** Purpose : converting 1-cat ice to 1 ice category 778 770 !!------------------------------------------------------------------- 779 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 780 REAL(wp), DIMENSION(:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 781 !!------------------------------------------------------------------- 782 zh_i(:) = zhti(:) 783 zh_s(:) = zhts(:) 784 za_i(:) = zati(:) 771 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 772 REAL(wp), DIMENSION(:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 773 REAL(wp), DIMENSION(:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds 774 REAL(wp), DIMENSION(:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds 775 !!------------------------------------------------------------------- 776 ! == thickness and concentration == ! 777 ph_i(:) = phti(:) 778 ph_s(:) = phts(:) 779 pa_i(:) = pati(:) 780 ! 781 ! == temperature and salinity and ponds == ! 782 pt_i (:) = ptmi (:) 783 pt_s (:) = ptms (:) 784 pt_su(:) = ptmsu(:) 785 ps_i (:) = psmi (:) 786 pa_ip(:) = patip(:) 787 ph_ip(:) = phtip(:) 788 785 789 END SUBROUTINE ice_var_itd_1c1c 786 790 787 SUBROUTINE ice_var_itd_Nc1c( zhti, zhts, zati, zh_i, zh_s, za_i ) 791 SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati , ph_i, ph_s, pa_i, & 792 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 788 793 !!------------------------------------------------------------------- 789 794 !! ** Purpose : converting N-cat ice to 1 ice category 790 795 !!------------------------------------------------------------------- 791 REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 792 REAL(wp), DIMENSION(:) , INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 793 !!------------------------------------------------------------------- 794 ! 795 za_i(:) = SUM( zati(:,:), dim=2 ) 796 ! 797 WHERE( za_i(:) /= 0._wp ) 798 zh_i(:) = SUM( zhti(:,:) * zati(:,:), dim=2 ) / za_i(:) 799 zh_s(:) = SUM( zhts(:,:) * zati(:,:), dim=2 ) / za_i(:) 800 ELSEWHERE 801 zh_i(:) = 0._wp 802 zh_s(:) = 0._wp 796 REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 797 REAL(wp), DIMENSION(:) , INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 798 REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds 799 REAL(wp), DIMENSION(:) , INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds 800 ! 801 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs 802 ! 803 INTEGER :: idim 804 !!------------------------------------------------------------------- 805 ! 806 idim = SIZE( phti, 1 ) 807 ! 808 ! == thickness and concentration == ! 809 ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim) ) 810 ! 811 pa_i(:) = SUM( pati(:,:), dim=2 ) 812 813 WHERE( ( pa_i(:) ) /= 0._wp ) ; z1_ai(:) = 1._wp / pa_i(:) 814 ELSEWHERE ; z1_ai(:) = 0._wp 803 815 END WHERE 816 817 ph_i(:) = SUM( phti(:,:) * pati(:,:), dim=2 ) * z1_ai(:) 818 ph_s(:) = SUM( phts(:,:) * pati(:,:), dim=2 ) * z1_ai(:) 819 ! 820 ! == temperature and salinity == ! 821 WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp ) ; z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) ) 822 ELSEWHERE ; z1_vi(:) = 0._wp 823 END WHERE 824 WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp ) ; z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) ) 825 ELSEWHERE ; z1_vs(:) = 0._wp 826 END WHERE 827 pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 828 pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 829 pt_su(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:) 830 ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 831 832 ! == ponds == ! 833 pa_ip(:) = SUM( patip(:,:), dim=2 ) 834 WHERE( pa_ip(:) /= 0._wp ) ; ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 835 ELSEWHERE ; ph_ip(:) = 0._wp 836 END WHERE 837 ! 838 DEALLOCATE( z1_ai, z1_vi, z1_vs ) 804 839 ! 805 840 END SUBROUTINE ice_var_itd_Nc1c 806 841 807 SUBROUTINE ice_var_itd_1cMc( zhti, zhts, zati, zh_i, zh_s, za_i ) 842 SUBROUTINE ice_var_itd_1cMc( phti, phts, pati , ph_i, ph_s, pa_i, & 843 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 808 844 !!------------------------------------------------------------------- 809 845 !! 810 846 !! ** Purpose : converting 1-cat ice to jpl ice categories 811 847 !! 812 !! ice thickness distribution follows a gaussian law 813 !! around the concentration of the most likely ice thickness 814 !! (similar as iceistate.F90) 815 !! 816 !! ** Method: Iterative procedure 817 !! 818 !! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian 819 !! 820 !! 2) Check whether the distribution conserves area and volume, positivity and 821 !! category boundaries 848 !! 849 !! ** Method: ice thickness distribution follows a gamma function from Abraham et al. (2015) 850 !! it has the property of conserving total concentration and volume 822 851 !! 823 !! 3) If not (input ice is too thin), the last category is empty and 824 !! the number of categories is reduced (jpl-1) 825 !! 826 !! 4) Iterate until ok (SUM(itest(:) = 4) 827 !! 828 !! ** Arguments : zhti: 1-cat ice thickness 829 !! zhts: 1-cat snow depth 830 !! zati: 1-cat ice concentration 852 !! 853 !! ** Arguments : phti: 1-cat ice thickness 854 !! phts: 1-cat snow depth 855 !! pati: 1-cat ice concentration 831 856 !! 832 857 !! ** Output : jpl-cat 833 858 !! 834 !! (Example of application: BDY forcings when input are cell averaged) 835 !!------------------------------------------------------------------- 836 INTEGER :: ji, jk, jl ! dummy loop indices 837 INTEGER :: idim, i_fill, jl0 838 REAL(wp) :: zarg, zV, zconv, zdh, zdv 839 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 840 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 841 INTEGER , DIMENSION(4) :: itest 842 !!------------------------------------------------------------------- 843 ! 844 ! ---------------------------------------- 845 ! distribution over the jpl ice categories 846 ! ---------------------------------------- 847 ! a gaussian distribution for ice concentration is used 848 ! then we check whether the distribution fullfills 849 ! volume and area conservation, positivity and ice categories bounds 850 idim = SIZE( zhti , 1 ) 851 zh_i(1:idim,1:jpl) = 0._wp 852 zh_s(1:idim,1:jpl) = 0._wp 853 za_i(1:idim,1:jpl) = 0._wp 854 ! 859 !! Abraham, C., Steiner, N., Monahan, A. and Michel, C., 2015. 860 !! Effects of subgrid‐scale snow thickness variability on radiative transfer in sea ice. 861 !! Journal of Geophysical Research: Oceans, 120(8), pp.5597-5614 862 !!------------------------------------------------------------------- 863 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 864 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 865 REAL(wp), DIMENSION(:) , INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds 866 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds 867 ! 868 REAL(wp), ALLOCATABLE, DIMENSION(:) :: zfra, z1_hti 869 INTEGER :: ji, jk, jl 870 INTEGER :: idim 871 REAL(wp) :: zv, zdh 872 !!------------------------------------------------------------------- 873 ! 874 idim = SIZE( phti , 1 ) 875 ! 876 ph_i(1:idim,1:jpl) = 0._wp 877 ph_s(1:idim,1:jpl) = 0._wp 878 pa_i(1:idim,1:jpl) = 0._wp 879 ! 880 ALLOCATE( z1_hti(idim) ) 881 WHERE( phti(:) /= 0._wp ) ; z1_hti(:) = 1._wp / phti(:) 882 ELSEWHERE ; z1_hti(:) = 0._wp 883 END WHERE 884 ! 885 ! == thickness and concentration == ! 886 ! for categories 1:jpl-1, integrate the gamma function from hi_max(jl-1) to hi_max(jl) 887 DO jl = 1, jpl-1 888 DO ji = 1, idim 889 ! 890 IF( phti(ji) > 0._wp ) THEN 891 ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl) 892 pa_i(ji,jl) = pati(ji) * z1_hti(ji) * ( ( phti(ji) + 2.*hi_max(jl-1) ) * EXP( -2.*hi_max(jl-1)*z1_hti(ji) ) & 893 & - ( phti(ji) + 2.*hi_max(jl ) ) * EXP( -2.*hi_max(jl )*z1_hti(ji) ) ) 894 ! 895 ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl) 896 zv = pati(ji) * z1_hti(ji) * ( ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jl-1) + 2.*hi_max(jl-1)*hi_max(jl-1) ) & 897 & * EXP( -2.*hi_max(jl-1)*z1_hti(ji) ) & 898 & - ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jl) + 2.*hi_max(jl)*hi_max(jl) ) & 899 & * EXP(-2.*hi_max(jl)*z1_hti(ji)) ) 900 ! thickness 901 IF( pa_i(ji,jl) > epsi06 ) THEN 902 ph_i(ji,jl) = zv / pa_i(ji,jl) 903 ELSE 904 ph_i(ji,jl) = 0. 905 pa_i(ji,jl) = 0. 906 ENDIF 907 ENDIF 908 ! 909 ENDDO 910 ENDDO 911 ! 912 ! for the last category (jpl), integrate the gamma function from hi_max(jpl-1) to infinity 855 913 DO ji = 1, idim 856 914 ! 857 IF( zhti(ji) > 0._wp ) THEN 858 ! 859 ! find which category (jl0) the input ice thickness falls into 860 jl0 = jpl 861 DO jl = 1, jpl 862 IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 863 jl0 = jl 864 CYCLE 865 ENDIF 866 END DO 867 ! 868 itest(:) = 0 869 i_fill = jpl + 1 !------------------------------------ 870 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 871 ! !------------------------------------ 872 i_fill = i_fill - 1 873 ! 874 zh_i(ji,1:jpl) = 0._wp 875 za_i(ji,1:jpl) = 0._wp 876 itest(:) = 0 877 ! 878 IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1 879 zh_i(ji,1) = zhti(ji) 880 za_i (ji,1) = zati (ji) 881 ELSE !-- case ice is thicker: fill categories >1 882 ! thickness 883 DO jl = 1, i_fill - 1 884 zh_i(ji,jl) = hi_mean(jl) 885 END DO 886 ! 887 ! concentration 888 za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl)) 889 DO jl = 1, i_fill - 1 890 IF ( jl /= jl0 ) THEN 891 zarg = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 892 za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2) 893 ENDIF 894 END DO 895 ! 896 ! last category 897 za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) ) 898 zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) ) 899 zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 900 ! 901 ! correction if concentration of upper cat is greater than lower cat 902 ! (it should be a gaussian around jl0 but sometimes it is not) 903 IF ( jl0 /= jpl ) THEN 904 DO jl = jpl, jl0+1, -1 905 IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN 906 zdv = zh_i(ji,jl) * za_i(ji,jl) 907 zh_i(ji,jl ) = 0._wp 908 za_i (ji,jl ) = 0._wp 909 za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 ) 910 END IF 911 END DO 912 ENDIF 913 ! 914 ENDIF 915 ! 916 ! Compatibility tests 917 zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 918 IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation 919 ! 920 zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) ) 921 IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation 922 ! 923 IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ? 924 ! 925 itest(4) = 1 926 DO jl = 1, i_fill 927 IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations 928 END DO 929 ! !---------------------------- 930 END DO ! end iteration on categories 931 ! !---------------------------- 915 IF( phti(ji) > 0._wp ) THEN 916 ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jpl-1) to infinity 917 pa_i(ji,jpl) = pati(ji) * z1_hti(ji) * ( phti(ji) + 2.*hi_max(jpl-1) ) * EXP( -2.*hi_max(jpl-1)*z1_hti(ji) ) 918 919 ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jpl-1) to infinity 920 zv = pati(ji) * z1_hti(ji) * ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jpl-1) + 2.*hi_max(jpl-1)*hi_max(jpl-1) ) & 921 & * EXP( -2.*hi_max(jpl-1)*z1_hti(ji) ) 922 ! thickness 923 IF( pa_i(ji,jpl) > epsi06 ) THEN 924 ph_i(ji,jpl) = zv / pa_i(ji,jpl) 925 else 926 ph_i(ji,jpl) = 0. 927 pa_i(ji,jpl) = 0. 928 ENDIF 932 929 ENDIF 933 END DO 934 935 ! Add Snow in each category where za_i is not 0 930 ! 931 ENDDO 932 ! 933 ! Add Snow in each category where pa_i is not 0 936 934 DO jl = 1, jpl 937 935 DO ji = 1, idim 938 IF( za_i(ji,jl) > 0._wp ) THEN939 zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji))936 IF( pa_i(ji,jl) > 0._wp ) THEN 937 ph_s(ji,jl) = ph_i(ji,jl) * phts(ji) * z1_hti(ji) 940 938 ! In case snow load is in excess that would lead to transformation from snow to ice 941 939 ! Then, transfer the snow excess into the ice (different from icethd_dh) 942 zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )940 zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rho0 ) * ph_i(ji,jl) ) * r1_rho0 ) 943 941 ! recompute h_i, h_s avoiding out of bounds values 944 zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )945 zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )942 ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh ) 943 ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos ) 946 944 ENDIF 947 945 END DO 948 946 END DO 949 947 ! 948 DEALLOCATE( z1_hti ) 949 ! 950 ! == temperature and salinity == ! 951 DO jl = 1, jpl 952 pt_i (:,jl) = ptmi (:) 953 pt_s (:,jl) = ptms (:) 954 pt_su(:,jl) = ptmsu(:) 955 ps_i (:,jl) = psmi (:) 956 ps_i (:,jl) = psmi (:) 957 END DO 958 ! 959 ! == ponds == ! 960 ALLOCATE( zfra(idim) ) 961 ! keep the same pond fraction atip/ati for each category 962 WHERE( pati(:) /= 0._wp ) ; zfra(:) = patip(:) / pati(:) 963 ELSEWHERE ; zfra(:) = 0._wp 964 END WHERE 965 DO jl = 1, jpl 966 pa_ip(:,jl) = zfra(:) * pa_i(:,jl) 967 END DO 968 ! keep the same v_ip/v_i ratio for each category 969 WHERE( ( phti(:) * pati(:) ) /= 0._wp ) ; zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) ) 970 ELSEWHERE ; zfra(:) = 0._wp 971 END WHERE 972 DO jl = 1, jpl 973 WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 974 ELSEWHERE ; ph_ip(:,jl) = 0._wp 975 END WHERE 976 END DO 977 DEALLOCATE( zfra ) 978 ! 950 979 END SUBROUTINE ice_var_itd_1cMc 951 980 952 SUBROUTINE ice_var_itd_NcMc( zhti, zhts, zati, zh_i, zh_s, za_i ) 981 SUBROUTINE ice_var_itd_NcMc( phti, phts, pati , ph_i, ph_s, pa_i, & 982 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 953 983 !!------------------------------------------------------------------- 954 984 !! … … 971 1001 !! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin) 972 1002 !! 973 !! ** Arguments : zhti: N-cat ice thickness974 !! zhts: N-cat snow depth975 !! zati: N-cat ice concentration1003 !! ** Arguments : phti: N-cat ice thickness 1004 !! phts: N-cat snow depth 1005 !! pati: N-cat ice concentration 976 1006 !! 977 1007 !! ** Output : jpl-cat … … 979 1009 !! (Example of application: BDY forcings when inputs have N-cat /= jpl) 980 1010 !!------------------------------------------------------------------- 981 INTEGER :: ji, jl, jl1, jl2 ! dummy loop indices 1011 REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 1012 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 1013 REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds 1014 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds 1015 ! 1016 INTEGER , ALLOCATABLE, DIMENSION(:,:) :: jlfil, jlfil2 1017 INTEGER , ALLOCATABLE, DIMENSION(:) :: jlmax, jlmin 1018 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs, ztmp, zfra 1019 ! 1020 REAL(wp), PARAMETER :: ztrans = 0.25_wp 1021 INTEGER :: ji, jl, jl1, jl2 982 1022 INTEGER :: idim, icat 983 REAL(wp), PARAMETER :: ztrans = 0.25_wp 984 REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 985 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 986 INTEGER , DIMENSION(:,:), ALLOCATABLE :: jlfil, jlfil2 987 INTEGER , DIMENSION(:) , ALLOCATABLE :: jlmax, jlmin 988 !!------------------------------------------------------------------- 989 ! 990 idim = SIZE( zhti, 1 ) 991 icat = SIZE( zhti, 2 ) 1023 !!------------------------------------------------------------------- 1024 ! 1025 idim = SIZE( phti, 1 ) 1026 icat = SIZE( phti, 2 ) 1027 ! 1028 ! == thickness and concentration == ! 992 1029 ! ! ---------------------- ! 993 1030 IF( icat == jpl ) THEN ! input cat = output cat ! 994 1031 ! ! ---------------------- ! 995 zh_i(:,:) = zhti(:,:) 996 zh_s(:,:) = zhts(:,:) 997 za_i(:,:) = zati(:,:) 1032 ph_i(:,:) = phti(:,:) 1033 ph_s(:,:) = phts(:,:) 1034 pa_i(:,:) = pati(:,:) 1035 ! 1036 ! == temperature and salinity and ponds == ! 1037 pt_i (:,:) = ptmi (:,:) 1038 pt_s (:,:) = ptms (:,:) 1039 pt_su(:,:) = ptmsu(:,:) 1040 ps_i (:,:) = psmi (:,:) 1041 pa_ip(:,:) = patip(:,:) 1042 ph_ip(:,:) = phtip(:,:) 998 1043 ! ! ---------------------- ! 999 1044 ELSEIF( icat == 1 ) THEN ! input cat = 1 ! 1000 1045 ! ! ---------------------- ! 1001 CALL ice_var_itd_1cMc( zhti(:,1), zhts(:,1), zati(:,1), zh_i(:,:), zh_s(:,:), za_i(:,:) ) 1046 CALL ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), & 1047 & ph_i(:,:), ph_s(:,:), pa_i (:,:), & 1048 & ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), & 1049 & pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:) ) 1002 1050 ! ! ---------------------- ! 1003 1051 ELSEIF( jpl == 1 ) THEN ! output cat = 1 ! 1004 1052 ! ! ---------------------- ! 1005 CALL ice_var_itd_Nc1c( zhti(:,:), zhts(:,:), zati(:,:), zh_i(:,1), zh_s(:,1), za_i(:,1) ) 1053 CALL ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), & 1054 & ph_i(:,1), ph_s(:,1), pa_i (:,1), & 1055 & ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), & 1056 & pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1) ) 1006 1057 ! ! ----------------------- ! 1007 1058 ELSE ! input cat /= output cat ! … … 1012 1063 1013 1064 ! --- initialize output fields to 0 --- ! 1014 zh_i(1:idim,1:jpl) = 0._wp1015 zh_s(1:idim,1:jpl) = 0._wp1016 za_i(1:idim,1:jpl) = 0._wp1065 ph_i(1:idim,1:jpl) = 0._wp 1066 ph_s(1:idim,1:jpl) = 0._wp 1067 pa_i(1:idim,1:jpl) = 0._wp 1017 1068 ! 1018 1069 ! --- fill the categories --- ! … … 1024 1075 DO jl2 = 1, icat 1025 1076 DO ji = 1, idim 1026 IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN1077 IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN 1027 1078 ! fill the right category 1028 zh_i(ji,jl1) = zhti(ji,jl2)1029 zh_s(ji,jl1) = zhts(ji,jl2)1030 za_i(ji,jl1) = zati(ji,jl2)1079 ph_i(ji,jl1) = phti(ji,jl2) 1080 ph_s(ji,jl1) = phts(ji,jl2) 1081 pa_i(ji,jl1) = pati(ji,jl2) 1031 1082 ! record categories that are filled 1032 1083 jlmax(ji) = MAX( jlmax(ji), jl1 ) … … 1045 1096 IF( jl1 > 1 ) THEN 1046 1097 ! fill the lower cat (jl1-1) 1047 za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)1048 zh_i(ji,jl1-1) = hi_mean(jl1-1)1098 pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1) 1099 ph_i(ji,jl1-1) = hi_mean(jl1-1) 1049 1100 ! remove from cat jl1 1050 za_i(ji,jl1 ) = ( 1._wp - ztrans ) * za_i(ji,jl1)1101 pa_i(ji,jl1 ) = ( 1._wp - ztrans ) * pa_i(ji,jl1) 1051 1102 ENDIF 1052 1103 IF( jl2 < jpl ) THEN 1053 1104 ! fill the upper cat (jl2+1) 1054 za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)1055 zh_i(ji,jl2+1) = hi_mean(jl2+1)1105 pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2) 1106 ph_i(ji,jl2+1) = hi_mean(jl2+1) 1056 1107 ! remove from cat jl2 1057 za_i(ji,jl2 ) = ( 1._wp - ztrans ) * za_i(ji,jl2)1108 pa_i(ji,jl2 ) = ( 1._wp - ztrans ) * pa_i(ji,jl2) 1058 1109 ENDIF 1059 1110 END DO … … 1065 1116 IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN 1066 1117 ! fill high 1067 za_i(ji,jl) = ztrans * za_i(ji,jl-1)1068 zh_i(ji,jl) = hi_mean(jl)1118 pa_i(ji,jl) = ztrans * pa_i(ji,jl-1) 1119 ph_i(ji,jl) = hi_mean(jl) 1069 1120 jlfil(ji,jl) = jl 1070 1121 ! remove low 1071 za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)1122 pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1) 1072 1123 ENDIF 1073 1124 END DO … … 1079 1130 IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN 1080 1131 ! fill low 1081 za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)1082 zh_i(ji,jl) = hi_mean(jl)1132 pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1) 1133 ph_i(ji,jl) = hi_mean(jl) 1083 1134 jlfil2(ji,jl) = jl 1084 1135 ! remove high 1085 za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)1136 pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1) 1086 1137 ENDIF 1087 1138 END DO … … 1090 1141 DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays 1091 1142 DEALLOCATE( jlmin, jlmax ) 1143 ! 1144 ! == temperature and salinity == ! 1145 ! 1146 ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) ) 1147 ! 1148 WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp ) ; z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 ) 1149 ELSEWHERE ; z1_ai(:) = 0._wp 1150 END WHERE 1151 WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp ) ; z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) 1152 ELSEWHERE ; z1_vi(:) = 0._wp 1153 END WHERE 1154 WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp ) ; z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) 1155 ELSEWHERE ; z1_vs(:) = 0._wp 1156 END WHERE 1157 ! 1158 ! fill all the categories with the same value 1159 ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 1160 DO jl = 1, jpl 1161 pt_i (:,jl) = ztmp(:) 1162 END DO 1163 ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 1164 DO jl = 1, jpl 1165 pt_s (:,jl) = ztmp(:) 1166 END DO 1167 ztmp(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:) 1168 DO jl = 1, jpl 1169 pt_su(:,jl) = ztmp(:) 1170 END DO 1171 ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 1172 DO jl = 1, jpl 1173 ps_i (:,jl) = ztmp(:) 1174 END DO 1175 ! 1176 DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp ) 1177 ! 1178 ! == ponds == ! 1179 ALLOCATE( zfra(idim) ) 1180 ! keep the same pond fraction atip/ati for each category 1181 WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp ) ; zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 ) 1182 ELSEWHERE ; zfra(:) = 0._wp 1183 END WHERE 1184 DO jl = 1, jpl 1185 pa_ip(:,jl) = zfra(:) * pa_i(:,jl) 1186 END DO 1187 ! keep the same v_ip/v_i ratio for each category 1188 WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp ) 1189 zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 ) 1190 ELSEWHERE 1191 zfra(:) = 0._wp 1192 END WHERE 1193 DO jl = 1, jpl 1194 WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 1195 ELSEWHERE ; ph_ip(:,jl) = 0._wp 1196 END WHERE 1197 END DO 1198 DEALLOCATE( zfra ) 1092 1199 ! 1093 1200 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
