- Timestamp:
- 2021-08-10T15:56:23+02:00 (3 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.4_FOAM_IAU_SI3_SIC/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.4_FOAM_IAU_SI3_SIC/src/ICE/icethd_do.F90
r14075 r15177 38 38 39 39 ! !!** namelist (namthd_do) ** 40 REAL(wp) :: rn_hinew! thickness for new ice formation (m)40 REAL(wp), PUBLIC :: rn_hinew ! thickness for new ice formation (m) 41 41 LOGICAL :: ln_frazil ! use of frazil ice collection as function of wind (T) or not (F) 42 42 REAL(wp) :: rn_maxfraz ! maximum portion of frazil ice collecting at the ice bottom -
NEMO/branches/UKMO/NEMO_4.0.4_FOAM_IAU_SI3_SIC/src/OCE/ASM/asminc.F90
r14075 r15177 36 36 USE diaobs , ONLY : calc_date ! Compute the calendar date on a given step 37 37 #if defined key_si3 38 USE ice , ONLY : hm_i, at_i, at_i_b 38 USE phycst ! physical constants 39 USE ice1D ! sea-ice: thermodynamics variables 40 USE icetab ! sea-ice: 2D <==> 1D 41 USE icethd_do 42 USE ice 43 USE icevar , ONLY : ice_var_zapsmall 39 44 #endif 40 45 ! … … 89 94 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ssh_bkg, ssh_bkginc ! Background sea surface height and its increment 90 95 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: seaice_bkginc ! Increment to the background sea ice conc 96 #if defined key_si3 97 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: a_i_bkginc ! Increment to the background sea ice conc categories 98 #endif 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 103 LOGICAL, PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: incr_newice !: mask .TRUE.=DA positive ice increment to open water 91 104 #if defined key_cice && defined key_asminc 92 105 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ndaice_da ! ice increment tendency into CICE … … 131 144 ! 132 145 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zhdiv ! 2D workspace 146 REAL(wp) :: zremaining_increment 147 148 #if defined key_si3 149 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z1_hti 150 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: hm_i_loc 151 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z1_at_i 152 REAL(wp) :: zopenwater_lim 153 INTEGER :: jl 154 #endif 133 155 !! 134 156 NAMELIST/nam_asminc/ ln_bkgwri, & 135 157 & ln_trainc, ln_dyninc, ln_sshinc, & 136 & ln_ asmdin, ln_asmiau, &158 & ln_seaiceinc, ln_asmdin, ln_asmiau, & 137 159 & nitbkg, nitdin, nitiaustr, nitiaufin, niaufn, & 138 160 & ln_salfix, salfixmin, nn_divdmp … … 144 166 ln_seaiceinc = .FALSE. 145 167 ln_temnofreeze = .FALSE. 168 na_iincr_split = 3 169 nh_iincr_min = 3 170 zhi_damin = 0.5_wp 171 zopenwater_lim = 1.0e-2_wp 146 172 147 173 REWIND( numnam_ref ) ! Namelist nam_asminc in reference namelist : Assimilation increment … … 173 199 WRITE(numout,*) ' Logical switch for ensuring that the sa > salfixmin ln_salfix = ', ln_salfix 174 200 WRITE(numout,*) ' Minimum salinity after applying the increments salfixmin = ', salfixmin 201 WRITE(numout,*) ' Ice splitting option na_iincr_split = ', na_iincr_split 202 WRITE(numout,*) ' Ice thickness for new ice from da option nh_iincr_min = ', nh_iincr_min 203 WRITE(numout,*) ' Minimum ice thickness for new ice from da zhi_damin = ', zhi_damin 204 WRITE(numout,*) ' Limit for open water detection for ice da zopenwater_lim = ', zopenwater_lim 175 205 ENDIF 176 206 … … 325 355 ALLOCATE( ssh_bkginc (jpi,jpj) ) ; ssh_bkginc (:,:) = 0._wp 326 356 ALLOCATE( seaice_bkginc(jpi,jpj) ) ; seaice_bkginc(:,:) = 0._wp 357 #if defined key_si3 358 ALLOCATE( a_i_bkginc (jpi,jpj,jpl) ) ; a_i_bkginc (:,:,:) = 0._wp 359 ALLOCATE( incr_newice(jpi,jpj,jpl) ) ; incr_newice(:,:,:) = .FALSE. 360 #endif 327 361 #if defined key_asminc 328 362 ALLOCATE( ssh_iau (jpi,jpj) ) ; ssh_iau (:,:) = 0._wp … … 398 432 ! to allow for differences in masks 399 433 WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0 434 435 IF (lwp) THEN 436 WRITE(numout,*) 437 WRITE(numout,*) 'asm_inc_init : Converting seaice_bkginc to a_i_bkginc' 438 WRITE(numout,*) '~~~~~~~~~~~~' 439 END IF 440 !single category increment for sea ice conc 441 !convert single category increment to multi category 442 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,*) '~~~~~~~~~~~~' 467 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(:,:) ) ) 487 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 400 570 ENDIF 401 571 ! … … 823 993 INTEGER, INTENT(in), OPTIONAL :: kindic ! flag for disabling the deallocation 824 994 ! 825 INTEGER :: it 995 INTEGER :: it, jl 826 996 REAL(wp) :: zincwgt ! IAU weight for current time step 827 997 #if defined key_si3 828 REAL(wp), DIMENSION(jpi,jpj) :: zofrld, zohicif, zseaicendg, zhicifinc 829 REAL(wp) :: zhicifmin = 0.5_wp ! ice minimum depth in metres 998 REAL(wp), DIMENSION(jpi,jpj,jpl) :: da_i ! change in ice concentration 999 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 metres 1001 REAL(wp) :: rn_hinew_save 1002 LOGICAL, SAVE :: initial_step=.TRUE. 1003 REAL(wp), DIMENSION(jpi,jpj) :: at_i_bkginc ! sum of conc increment over categories 830 1004 #endif 831 1005 !!---------------------------------------------------------------------- … … 850 1024 ! 851 1025 #if defined key_si3 852 zofrld (:,:) = 1._wp - at_i(:,:) 853 zohicif(:,:) = hm_i(:,:) 854 ! 855 at_i (:,:) = 1. - MIN( MAX( 1.-at_i (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 856 at_i_b(:,:) = 1. - MIN( MAX( 1.-at_i_b(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 857 fr_i(:,:) = at_i(:,:) ! adjust ice fraction 858 ! 859 zseaicendg(:,:) = zofrld(:,:) - (1. - at_i(:,:)) ! find out actual sea ice nudge applied 860 ! 861 ! Nudge sea ice depth to bring it up to a required minimum depth 862 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 863 zhicifinc(:,:) = (zhicifmin - hm_i(:,:)) * zincwgt 864 ELSEWHERE 865 zhicifinc(:,:) = 0.0_wp 866 END WHERE 867 ! 868 ! nudge ice depth 869 hm_i (:,:) = hm_i (:,:) + zhicifinc(:,:) 1026 IF ( zhi_damin > hi_max(1) - epsi10 ) THEN 1027 zhi_damin = hi_max(1) - epsi10 1028 ENDIF 1029 1030 where( a_i(:,:,:) > epsi10 ) 1031 z1_a_i(:,:,:) = 1.0_wp/a_i(:,:,:) 1032 elsewhere 1033 z1_a_i(:,:,:) = 0.0_wp 1034 end where 1035 1036 ! add concentration increments and bound above zero 1037 where (.not. incr_newice) 1038 a_i(:,:,:) = MAX(a_i(:,:,:) + a_i_bkginc(:,:,:) * zincwgt, 0.0_wp) 1039 end where 1040 1041 ! ensure total sum of categories doesn't exceed rn_amax_2d 1042 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 1043 DO jl = 1, jpl 1044 WHERE( at_i(:,:) > rn_amax_2d(:,:) ) a_i(:,:,jl) = a_i(:,:,jl) * rn_amax_2d(:,:) / at_i(:,:) 1045 END DO 1046 1047 ! compute change in ice concentration (new / old) 1048 where (incr_newice) 1049 da_i(:,:,:) = 1.0_wp 1050 elsewhere 1051 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 1141 1142 ! make volume consistent with concentration and thickness 1143 ! note where conc is zero this will lead to zero volume 1144 ! and soon zero thickness after the call to ice_var_zapsmall 1145 v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:) 1146 1147 ! as ice concentration has changed, maintain equivalent values for prognostic variables 1148 v_s (:,:,:) = v_s (:,:,:) * da_i(:,:,:) 1149 sv_i (:,:,:) = sv_i(:,:,:) * da_i(:,:,:) 1150 a_ip (:,:,:) = a_ip(:,:,:) * da_i(:,:,:) 1151 v_ip (:,:,:) = v_ip(:,:,:) * da_i(:,:,:) 1152 do jl = 1, nlay_i 1153 e_i(:,:,jl,:) = e_i(:,:,jl,:) * da_i(:,:,:) 1154 end do 1155 do jl = 1, nlay_s 1156 e_s(:,:,jl,:) = e_s(:,:,jl,:) * da_i(:,:,:) 1157 end do 1158 1159 call ice_var_zapsmall() 1160 870 1161 ! 871 1162 ! seaice salinity balancing (to add) … … 895 1186 ! 896 1187 neuler = 0 ! Force Euler forward step 897 ! 898 ! Sea-ice : SI3 case 899 ! 900 #if defined key_si3 901 zofrld (:,:) = 1._wp - at_i(:,:) 902 zohicif(:,:) = hm_i(:,:) 903 ! 904 ! Initialize the now fields the background + increment 905 at_i(:,:) = 1. - MIN( MAX( 1.-at_i(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp) 906 at_i_b(:,:) = at_i(:,:) 907 fr_i(:,:) = at_i(:,:) ! adjust ice fraction 908 ! 909 zseaicendg(:,:) = zofrld(:,:) - (1. - at_i(:,:)) ! find out actual sea ice nudge applied 910 ! 911 ! Nudge sea ice depth to bring it up to a required minimum depth 912 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 913 zhicifinc(:,:) = zhicifmin - hm_i(:,:) 914 ELSEWHERE 915 zhicifinc(:,:) = 0.0_wp 916 END WHERE 917 ! 918 ! nudge ice depth 919 hm_i (:,:) = hm_i (:,:) + zhicifinc(:,:) 920 ! 921 ! seaice salinity balancing (to add) 922 #endif 1188 1189 IF(lwp) THEN 1190 WRITE(numout,*) 1191 WRITE(numout,*) 'seaice_asm_inc : sea ice direct initialization at time step = ', kt 1192 WRITE(numout,*) '~~~~~~~~~~~~' 1193 ENDIF 923 1194 ! 924 1195 #if defined key_cice && defined key_asminc … … 937 1208 ! 938 1209 ENDIF 939 940 !#if defined defined key_si3 || defined key_cice941 !942 ! IF (ln_seaicebal ) THEN943 ! !! balancing salinity increments944 ! !! simple case from limflx.F90 (doesn't include a mass flux)945 ! !! assumption is that as ice concentration is reduced or increased946 ! !! the snow and ice depths remain constant947 ! !! note that snow is being created where ice concentration is being increased948 ! !! - could be more sophisticated and949 ! !! not do this (but would need to alter h_snow)950 !951 ! usave(:,:,:)=sb(:,:,:) ! use array as a temporary store952 !953 ! DO jj = 1, jpj954 ! DO ji = 1, jpi955 ! ! calculate change in ice and snow mass per unit area956 ! ! positive values imply adding salt to the ocean (results from ice formation)957 ! ! fwf : ice formation and melting958 !959 ! zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt960 !961 ! ! change salinity down to mixed layer depth962 ! mld=hmld_kara(ji,jj)963 !964 ! ! prevent small mld965 ! ! less than 10m can cause salinity instability966 ! IF (mld < 10) mld=10967 !968 ! ! set to bottom of a level969 ! DO jk = jpk-1, 2, -1970 ! IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN971 ! mld=gdepw(ji,jj,jk+1)972 ! jkmax=jk973 ! ENDIF974 ! ENDDO975 !976 ! ! avoid applying salinity balancing in shallow water or on land977 ! !978 !979 ! ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)980 !981 ! dsal_ocn=0.0_wp982 ! sal_thresh=5.0_wp ! minimum salinity threshold for salinity balancing983 !984 ! if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) &985 ! dsal_ocn = zfons / (rhop(ji,jj,1) * mld)986 !987 ! ! put increments in for levels in the mixed layer988 ! ! but prevent salinity below a threshold value989 !990 ! DO jk = 1, jkmax991 !992 ! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN993 ! sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn994 ! sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn995 ! ENDIF996 !997 ! ENDDO998 !999 ! ! ! salt exchanges at the ice/ocean interface1000 ! ! zpmess = zfons / rdt_ice ! rdt_ice is ice timestep1001 ! !1002 ! !! Adjust fsalt. A +ve fsalt means adding salt to ocean1003 ! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt1004 ! !!1005 ! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d)1006 ! !! ! E-P (kg m-2 s-2)1007 ! ! emp(ji,jj) = emp(ji,jj) + zpmess ! E-P (kg m-2 s-2)1008 ! ENDDO !ji1009 ! ENDDO !jj!1010 !1011 ! ENDIF !ln_seaicebal1012 !1013 !#endif1014 1210 ! 1015 1211 ENDIF … … 1017 1213 END SUBROUTINE seaice_asm_inc 1018 1214 1215 SUBROUTINE seaice_asm_qlead(ht_i_new, at_i_new, qlead_new) 1216 !!---------------------------------------------------------------------- 1217 !! *** ROUTINE seaice_asm_qlead *** 1218 !! 1219 !! ** Purpose : Compute the value of qlead that will correspond to new ice of 1220 !! thickness ht_i_new concentration of at_i_new 1221 !! 1222 !! ** Method : Based on symbolic inversion of the icethd_do routine 1223 !! 1224 !! ** Action : return qlead_new 1225 !!---------------------------------------------------------------------- 1226 REAL(wp), DIMENSION(:,:), INTENT(in) :: ht_i_new ! total thickness of new ice 1227 REAL(wp), DIMENSION(:,:), INTENT(in) :: at_i_new ! total concentration of new ice 1228 REAL(wp), DIMENSION(:,:), INTENT(inout) :: qlead_new ! output flux? value 1229 1230 INTEGER :: jj, ji 1231 1232 REAL(wp), DIMENSION(jpij) :: at_i_new_1d ! 1d version of at_i_new 1233 REAL(wp), DIMENSION(jpij) :: zs_newice ! salinity of accreted ice 1234 REAL(wp), DIMENSION(jpij) :: zh_newice ! thickness of accreted ice 1235 !!---------------------------------------------------------------------- 1236 1237 ! Identify grid points where new ice forms 1238 npti = 0 ; nptidx(:) = 0 1239 DO jj = 1, jpj 1240 DO ji = 1, jpi 1241 IF ( ht_i_new(ji,jj) > 0._wp ) THEN 1242 npti = npti + 1 1243 nptidx( npti ) = (jj - 1) * jpi + ji 1244 ENDIF 1245 END DO 1246 END DO 1247 1248 ! Move from 2-D to 1-D vectors 1249 IF ( npti > 0 ) THEN 1250 CALL tab_2d_1d( npti, nptidx(1:npti), at_i_new_1d(1:npti), at_i_new ) 1251 CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , ht_i_new ) 1252 CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d (1:npti) , sss_m ) 1253 CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d (1:npti) , qlead_new ) 1254 CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d (1:npti) , t_bo ) 1255 1256 1257 ! --- Salinity of new ice --- ! 1258 SELECT CASE ( nn_icesal ) 1259 CASE ( 1 ) ! Sice = constant 1260 zs_newice(1:npti) = rn_icesal 1261 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 1262 DO ji = 1, npti 1263 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) ) 1264 END DO 1265 CASE ( 3 ) ! Sice = F(z) [multiyear ice] 1266 zs_newice(1:npti) = 2.3 1267 END SELECT 1268 1269 1270 DO ji = 1, npti 1271 qlead_1d(ji) = at_i_new_1d(ji)*zh_newice(ji)*(-r1_rhoi*rLfus*rTmlt*rhoi*zs_newice(ji) + & 1272 & r1_rhoi*rhoi*(-rLfus - rTmlt*rcp*zs_newice(ji) + rTmlt*rcpi*zs_newice(ji) - & 1273 & rcpi*rt0 + rcpi*t_bo_1d(ji))* & 1274 & min(-epsi10, -rt0 + t_bo_1d(ji)) + & 1275 & rcp*(rt0 - t_bo_1d(ji))*min(-epsi10, -rt0 + t_bo_1d(ji))) / & 1276 & (r1_rhoi*min(-epsi10, -rt0 + t_bo_1d(ji))) 1277 END DO 1278 CALL tab_1d_2d( npti, nptidx(1:npti), qlead_1d(1:npti), qlead_new ) 1279 ENDIF ! npti > 0 1280 END SUBROUTINE seaice_asm_qlead 1281 1019 1282 !!====================================================================== 1020 1283 END MODULE asminc
Note: See TracChangeset
for help on using the changeset viewer.