Changeset 15241 for NEMO/branches
- Timestamp:
- 2021-09-10T10:44:30+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.4_FOAM_IAU_SI3_SIC/src/OCE/ASM/asminc.F90
r15177 r15241 98 98 #endif 99 99 REAL(wp) :: zhi_damin = 0.5_wp !: ice thickness for new sea ice from da increment 100 INTEGER :: na_iincr_split = 0 !: methodology for splitting a_i increment 101 INTEGER :: nh_iincr_min = 0 !: methodology for setting minimum ice thickness with DA 102 INTEGER :: nhi_damin = 1 !: thickness category corresponding to zhi_damin 100 INTEGER :: nhi_damin !: thickness category corresponding to zhi_damin 103 101 LOGICAL, PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: incr_newice !: mask .TRUE.=DA positive ice increment to open water 104 102 #if defined key_cice && defined key_asminc … … 147 145 148 146 #if defined key_si3 149 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z1_hti150 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: hm_i_loc151 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z1_at_i152 147 REAL(wp) :: zopenwater_lim 153 148 INTEGER :: jl … … 166 161 ln_seaiceinc = .FALSE. 167 162 ln_temnofreeze = .FALSE. 168 na_iincr_split = 3169 nh_iincr_min = 3170 163 zhi_damin = 0.5_wp 171 164 zopenwater_lim = 1.0e-2_wp … … 199 192 WRITE(numout,*) ' Logical switch for ensuring that the sa > salfixmin ln_salfix = ', ln_salfix 200 193 WRITE(numout,*) ' Minimum salinity after applying the increments salfixmin = ', salfixmin 201 WRITE(numout,*) ' Ice splitting option na_iincr_split = ', na_iincr_split202 WRITE(numout,*) ' Ice thickness for new ice from da option nh_iincr_min = ', nh_iincr_min203 194 WRITE(numout,*) ' Minimum ice thickness for new ice from da zhi_damin = ', zhi_damin 204 195 WRITE(numout,*) ' Limit for open water detection for ice da zopenwater_lim = ', zopenwater_lim … … 441 432 !convert single category increment to multi category 442 433 a_i_bkginc = 0.0_wp 443 444 ! hm_i seems to be inaccurate in areas of low ice conc 445 ! so compute our own estimate less prone to numerical issues 446 447 !ensure total volume and conc are up-to-date 448 vt_i = sum(v_i(:,:,:), dim=3) 449 at_i = sum(a_i(:,:,:), dim=3) 450 451 ALLOCATE( z1_hti , MOLD=vt_i ) 452 ALLOCATE( hm_i_loc, MOLD=vt_i ) 453 ALLOCATE( z1_at_i , MOLD=vt_i ) 454 455 nhi_damin=1 456 IF ( zhi_damin > hi_max(1) - epsi10 ) THEN 457 zhi_damin = hi_max(1) - epsi10 458 ENDIF 459 460 select case(na_iincr_split) 461 462 case( 0 ) ! gamma function 463 IF (lwp) THEN 464 WRITE(numout,*) 465 WRITE(numout,*) 'asm_inc_init : Converting seaice_bkginc to a_i_bkginc using gamma function' 466 WRITE(numout,*) '~~~~~~~~~~~~' 434 at_i = SUM(a_i(:,:,:), DIM=3) 435 436 ! Calculate which category corresponds to zhi_damin 437 ! find which category to fill 438 DO jl = 1, jpl 439 IF( zhi_damin > hi_max(jl-1) .AND. zhi_damin <= hi_max(jl) ) THEN 440 nhi_damin = jl 467 441 END IF 468 469 ! compute the inverse of total conc (cautiously) 470 where( at_i(:,:) > 0.01_wp) 471 z1_at_i(:,:) = 1.0_wp/at_i(:,:) 472 elsewhere 473 z1_at_i(:,:) = 0.0_wp 474 end where 475 476 ! compute mean thickness of background- call it hm_i_loc 477 hm_i_loc(:,:) = vt_i(:,:)*z1_at_i(:,:) 478 hm_i_loc(:,:) = MAX(hm_i_loc(:,:), 0.2_wp) 479 z1_hti(:,:) = 1._wp / hm_i_loc(:,:) 480 481 ! == thickness and concentration == ! 482 ! for categories 1:jpl-1, integrate the gamma function from hi_max(jl-1) to hi_max(jl) 483 DO jl = 1, jpl-1 484 ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl) 485 a_i_bkginc(:,:,jl) = seaice_bkginc(:,:) * z1_hti(:,:) * ( ( hm_i_loc(:,:) + 2.*hi_max(jl-1) ) * EXP( -2.*hi_max(jl-1)*z1_hti(:,:) ) & 486 & - ( hm_i_loc(:,:) + 2.*hi_max(jl ) ) * EXP( -2.*hi_max(jl )*z1_hti(:,:) ) ) 442 END DO 443 444 IF (lwp) THEN 445 WRITE(numout,*) 446 WRITE(numout,*) 'asm_inc_init : Converting seaice_bkginc to a_i_bkginc using Peterson splitting' 447 WRITE(numout,*) '~~~~~~~~~~~~' 448 END IF 449 DO jj = 1, jpj 450 DO ji = 1, jpi 451 IF ( seaice_bkginc(ji,jj) > 0.0_wp) THEN 452 !Positive ice concentration increments are always 453 !added to the thinnest category of ice 454 a_i_bkginc(ji,jj,nhi_damin) = seaice_bkginc(ji,jj) 455 ELSE 456 !negative increments are first removed from the thinnest 457 !available category until it reaches zero concentration 458 !and then progressively removed from thicker categories 459 zremaining_increment = seaice_bkginc(ji,jj) 460 DO jl = 1, jpl 461 ! assign as much increment as possible to current category 462 a_i_bkginc(ji,jj,jl) = -min(a_i(ji,jj,jl), -zremaining_increment) 463 ! update remaining amount of increment 464 zremaining_increment = zremaining_increment - a_i_bkginc(ji,jj,jl) 465 END DO 466 END IF 487 467 END DO 488 ! for the last category (jpl), integrate the gamma function from hi_max(jpl-1) to infinity 489 where ( hm_i_loc(:,:) > 0._wp ) 490 ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jpl-1) to infinity 491 a_i_bkginc(:,:,jpl) = seaice_bkginc(:,:) * z1_hti(:,:) * ( hm_i_loc(:,:) + 2.*hi_max(jpl-1) ) * EXP( -2.*hi_max(jpl-1)*z1_hti(:,:) ) 492 END where 493 494 case (1) ! maximum a_i location 495 IF (lwp) THEN 496 WRITE(numout,*) 497 WRITE(numout,*) 'asm_inc_init : Converting seaice_bkginc to a_i_bkginc using max a_i location' 498 WRITE(numout,*) '~~~~~~~~~~~~' 499 END IF 500 DO jj = 1, jpj 501 DO ji = 1, jpi 502 if (at_i(ji,jj) > 0.01_wp) then 503 jl = maxloc(pack(a_i(ji,jj,:), .true.), dim=1) 504 else 505 jl = nhi_damin 506 endif 507 a_i_bkginc(ji,jj,jl) = seaice_bkginc(ji,jj) 508 end do 509 end do 510 case (2) ! background profile 511 IF (lwp) THEN 512 WRITE(numout,*) 513 WRITE(numout,*) 'asm_inc_init : Converting seaice_bkginc to a_i_bkginc using background profile' 514 WRITE(numout,*) '~~~~~~~~~~~~' 515 END IF 516 DO jj = 1, jpj 517 DO ji = 1, jpi 518 if (at_i(ji,jj) > 0.01_wp) then 519 do jl = 1, jpl 520 a_i_bkginc(ji,jj,jl) = a_i(ji,jj,jl) * seaice_bkginc(ji,jj) / at_i(ji,jj) 521 end do 522 else 523 a_i_bkginc(ji,jj,nhi_damin) = seaice_bkginc(ji,jj) 524 end if 525 end do 526 end do 527 case (3) ! UK MO scheme based on Peterson et al. 2015 528 ! https://doi.org/10.1007/s00382-014-2190-9 529 ! but implemented by PB at ECMWF 530 IF (lwp) THEN 531 WRITE(numout,*) 532 WRITE(numout,*) 'asm_inc_init : Converting seaice_bkginc to a_i_bkginc using Peterson splitting' 533 WRITE(numout,*) '~~~~~~~~~~~~' 534 END IF 535 DO jj = 1, jpj 536 DO ji = 1, jpi 537 IF ( seaice_bkginc(ji,jj) > 0.0_wp) THEN 538 !Positive ice concentration increments are always 539 !added to the thinnest category of ice 540 a_i_bkginc(ji,jj,nhi_damin) = seaice_bkginc(ji,jj) 541 ELSE 542 !negative increments are first removed from the thinnest 543 !available category until it reaches zero concentration 544 !and then progressively removed from thicker categories 545 zremaining_increment = seaice_bkginc(ji,jj) 546 DO jl = 1, jpl 547 ! assign as much increment as possible to current category 548 a_i_bkginc(ji,jj,jl) = -min(a_i(ji,jj,jl), -zremaining_increment) 549 ! update remaining amount of increment 550 zremaining_increment = zremaining_increment - a_i_bkginc(ji,jj,jl) 551 END DO 552 END IF 553 END DO 554 END DO 555 case default 556 CALL ctl_stop('Bad choice of na_iincr_split') 557 end select 558 559 do jl = 1,jpl 560 where (at_i(:,:) < zopenwater_lim .and. seaice_bkginc(:,:) > 0.0_wp) 561 incr_newice(:,:,jl) = .true. 562 end where 563 end do 564 565 566 deallocate(z1_hti) 567 deallocate(z1_at_i) 568 deallocate(hm_i_loc) 569 468 END DO 469 DO jl = 1,jpl 470 WHERE (at_i(:,:) < zopenwater_lim .AND. seaice_bkginc(:,:) > 0.0_wp) 471 incr_newice(:,:,jl) = .TRUE. 472 END WHERE 473 END DO 570 474 ENDIF 571 475 ! … … 998 902 REAL(wp), DIMENSION(jpi,jpj,jpl) :: da_i ! change in ice concentration 999 903 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i ! inverse of old ice concentration 1000 REAL(wp) :: zhicifmin = 0.5_wp ! ice minimum depth in metres1001 904 REAL(wp) :: rn_hinew_save 1002 905 LOGICAL, SAVE :: initial_step=.TRUE. … … 1024 927 ! 1025 928 #if defined key_si3 1026 IF ( zhi_damin > hi_max(1) - epsi10 ) THEN 1027 zhi_damin = hi_max(1) - epsi10 1028 ENDIF 1029 1030 where( a_i(:,:,:) > epsi10 ) 929 WHERE( a_i(:,:,:) > epsi10 ) 1031 930 z1_a_i(:,:,:) = 1.0_wp/a_i(:,:,:) 1032 elsewhere931 ELSEWHERE 1033 932 z1_a_i(:,:,:) = 0.0_wp 1034 end where 1035 1036 ! add concentration increments and bound above zero 1037 where (.not. incr_newice) 933 END WHERE 934 935 WRITE(numout,*) "zhi_damin = ",zhi_damin 936 ! add positive concentration increments with zhi_damin to regions where ice 937 ! is already present and bound them to 1 938 ! ice volume is added based on the sea ice conc increment and assigned thickness 939 WHERE ( .NOT. incr_newice .AND. a_i_bkginc(:,:,:) > 0.0_wp ) 940 a_i(:,:,:) = a_i(:,:,:) + MIN(1.0_wp - a_i(:,:,:), a_i_bkginc(:,:,:) * zincwgt) 941 v_i(:,:,:) = v_i(:,:,:) + MIN(1.0_wp - a_i(:,:,:), a_i_bkginc(:,:,:) * zincwgt) * zhi_damin 942 END WHERE 943 944 ! add negative concentration increments to regions where ice 945 ! is already present and bound them to 0 946 ! in this case ice volume is changed based on the current thickness 947 WHERE ( .NOT. incr_newice .AND. a_i_bkginc(:,:,:) < 0.0_wp ) 1038 948 a_i(:,:,:) = MAX(a_i(:,:,:) + a_i_bkginc(:,:,:) * zincwgt, 0.0_wp) 1039 end where 949 v_i(:,:,:) = a_i(:,:,:) * h_i(:,:,:) 950 END WHERE 1040 951 1041 ! ensure total sum of categories doesn't exceed rn_amax_2d1042 at_i(:,:) = SUM( a_i(:,:,:), dim=3 )1043 DO jl = 1, jpl1044 WHERE( at_i(:,:) > rn_amax_2d(:,:) ) a_i(:,:,jl) = a_i(:,:,jl) * rn_amax_2d(:,:) / at_i(:,:)1045 END DO1046 1047 952 ! compute change in ice concentration (new / old) 1048 where (incr_newice)953 WHERE ( incr_newice ) 1049 954 da_i(:,:,:) = 1.0_wp 1050 elsewhere955 ELSEWHERE 1051 956 da_i(:,:,:) = a_i(:,:,:) * z1_a_i(:,:,:) 1052 end where 1053 1054 ! add thickness increments here 1055 1056 select case(nh_iincr_min) 1057 case (0) 1058 ! ensure minimum thickness where concentration is more than zero 1059 ! and maybe snow thickness? 1060 where ( a_i (:,:,:) > 0.001_wp ) 1061 h_i (:,:,:) = max(h_i (:,:,:), 0.01_wp) 1062 !v_s (:,:,:) = max(v_s (:,:,:), 0.01_wp) ! not a clue if this is a sensible value 1063 end where 1064 case (1) 1065 ! set minimum thickness only when increment has 1066 ! added ice to previously ice free points 1067 where ( z1_a_i(:,:,:) == 0.0_wp .and. & !old ice was zero 1068 & a_i(:,:,:) > epsi10 ) !new ice added 1069 h_i (:,:,:) = zhi_damin 1070 end where 1071 case (2) 1072 ! set minimum thickness only when increment has 1073 ! added ice to previously ice free points 1074 where ( z1_a_i(:,:,:) == 0.0_wp .and. & !old ice was zero 1075 & a_i(:,:,:) > epsi10 ) !new ice added 1076 h_i (:,:,:) = zhi_damin*hi_max(1) 1077 end where 1078 case (3) 1079 ! compute heat balance that adds specified ice thickness 1080 ! to open water 1081 if (initial_step) then 1082 IF(lwp) THEN 1083 WRITE(numout,*) 1084 WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' compute heat balance qlead' 1085 WRITE(numout,*) '~~~~~~~~~~~~' 1086 ENDIF 1087 1088 ! compute qlead to ensure thickness of zhi_damin 1089 ! for the new ice concentration values 1090 where (any(incr_newice, dim=3)) 1091 ht_i_new(:,:) = 0.1_wp 1092 elsewhere 1093 ht_i_new(:,:) = 0.0_wp 1094 end where 1095 1096 ! ensure all categories of new ice are zero 1097 where ( incr_newice ) 1098 v_i (:,:,:) = 0.0_wp 1099 v_s (:,:,:) = 0.0_wp 1100 sv_i(:,:,:) = 0.0_wp 1101 a_ip(:,:,:) = 0.0_wp 1102 v_ip(:,:,:) = 0.0_wp 1103 end where 1104 do jl = 1, nlay_i 1105 where ( incr_newice ) 1106 e_i(:,:,jl,:) = 0.0_wp 1107 end where 1108 end do 1109 do jl = 1, nlay_s 1110 where ( incr_newice ) 1111 e_s(:,:,jl,:) = 0.0_wp 1112 end where 1113 end do 1114 1115 qlead = 0.0_wp 1116 at_i_bkginc (:,:) = SUM( a_i_bkginc(:,:,:), dim=3 ) 1117 1118 call seaice_asm_qlead(ht_i_new, at_i_bkginc, qlead) 1119 1120 ! Call ice_thd_do to create the new ice from open water 1121 ! Note this adds the entire concentration increment too 1122 rn_hinew_save = rn_hinew 1123 rn_hinew = zhi_damin 1124 call ice_thd_do 1125 rn_hinew = rn_hinew_save 1126 1127 ! compute thickness now 1128 where( a_i(:,:,:) > epsi10 ) 1129 z1_a_i(:,:,:) = 1.0_wp/a_i(:,:,:) 1130 elsewhere 1131 z1_a_i(:,:,:) = 0.0_wp 1132 end where 1133 h_i(:,:,:) = v_i(:,:,:) * z1_a_i(:,:,:) 1134 1135 initial_step = .false. 1136 end if 1137 case default 1138 CALL ctl_stop('Bad choice of nh_iincr_min') 1139 end select 1140 957 END WHERE 958 959 ! compute heat balance that adds specified ice thickness 960 ! to open water 961 IF (initial_step) THEN 962 IF(lwp) THEN 963 WRITE(numout,*) 964 WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' compute heat balance qlead' 965 WRITE(numout,*) '~~~~~~~~~~~~' 966 ENDIF 967 968 ! compute qlead to ensure thickness of zhi_damin 969 ! for the new ice concentration values 970 WHERE (ANY(incr_newice, DIM=3)) 971 ht_i_new(:,:) = zhi_damin 972 ELSEWHERE 973 ht_i_new(:,:) = 0.0_wp 974 ENDWHERE 975 976 ! ensure all categories of new ice are zero 977 WHERE ( incr_newice ) 978 v_i (:,:,:) = 0.0_wp 979 v_s (:,:,:) = 0.0_wp 980 sv_i(:,:,:) = 0.0_wp 981 a_ip(:,:,:) = 0.0_wp 982 v_ip(:,:,:) = 0.0_wp 983 END WHERE 984 DO jl = 1, nlay_i 985 WHERE ( incr_newice ) 986 e_i(:,:,jl,:) = 0.0_wp 987 END WHERE 988 END DO 989 DO jl = 1, nlay_s 990 WHERE ( incr_newice ) 991 e_s(:,:,jl,:) = 0.0_wp 992 END WHERE 993 END DO 994 995 qlead = 0.0_wp 996 at_i_bkginc(:,:) = SUM( a_i_bkginc(:,:,:)*zincwgt, DIM=3 ) 997 998 call seaice_asm_qlead(ht_i_new, at_i_bkginc, qlead) 999 1000 ! Call ice_thd_do to create the new ice from open water 1001 rn_hinew_save = rn_hinew 1002 rn_hinew = zhi_damin 1003 call ice_thd_do 1004 rn_hinew = rn_hinew_save 1005 1006 ! compute thickness now 1007 WHERE( a_i(:,:,:) > epsi10 ) 1008 z1_a_i(:,:,:) = 1.0_wp/a_i(:,:,:) 1009 ELSEWHERE 1010 z1_a_i(:,:,:) = 0.0_wp 1011 END WHERE 1012 h_i(:,:,:) = v_i(:,:,:) * z1_a_i(:,:,:) 1013 1014 incr_newice(:,:,:) = .FALSE. 1015 initial_step = .FALSE. 1016 END IF 1141 1017 1142 1018 ! make volume consistent with concentration and thickness … … 1150 1026 a_ip (:,:,:) = a_ip(:,:,:) * da_i(:,:,:) 1151 1027 v_ip (:,:,:) = v_ip(:,:,:) * da_i(:,:,:) 1152 dojl = 1, nlay_i1028 DO jl = 1, nlay_i 1153 1029 e_i(:,:,jl,:) = e_i(:,:,jl,:) * da_i(:,:,:) 1154 end do1155 dojl = 1, nlay_s1030 END DO 1031 DO jl = 1, nlay_s 1156 1032 e_s(:,:,jl,:) = e_s(:,:,jl,:) * da_i(:,:,:) 1157 end do1033 END DO 1158 1034 1159 1035 call ice_var_zapsmall()
Note: See TracChangeset
for help on using the changeset viewer.
