New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 15177 – NEMO

Changeset 15177


Ignore:
Timestamp:
2021-08-10T15:56:23+02:00 (3 years ago)
Author:
dcarneir
Message:

Options for SIC IAU with SI3

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  
    3838 
    3939   !                          !!** 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) 
    4141   LOGICAL  ::   ln_frazil     ! use of frazil ice collection as function of wind (T) or not (F) 
    4242   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  
    3636   USE diaobs   , ONLY : calc_date     ! Compute the calendar date on a given step 
    3737#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 
    3944#endif 
    4045   ! 
    … …  
    8994   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment 
    9095   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 
    91104#if defined key_cice && defined key_asminc 
    92105   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ndaice_da             ! ice increment tendency into CICE 
    … …  
    131144      ! 
    132145      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 
    133155      !! 
    134156      NAMELIST/nam_asminc/ ln_bkgwri,                                      & 
    135157         &                 ln_trainc, ln_dyninc, ln_sshinc,                & 
    136          &                 ln_asmdin, ln_asmiau,                           & 
     158         &                 ln_seaiceinc, ln_asmdin, ln_asmiau,                           & 
    137159         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   & 
    138160         &                 ln_salfix, salfixmin, nn_divdmp 
    … …  
    144166      ln_seaiceinc   = .FALSE. 
    145167      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 
    146172 
    147173      REWIND( numnam_ref )              ! Namelist nam_asminc in reference namelist : Assimilation increment 
    … …  
    173199         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix 
    174200         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 
    175205      ENDIF 
    176206 
    … …  
    325355      ALLOCATE( ssh_bkginc   (jpi,jpj)     )   ;   ssh_bkginc   (:,:)   = 0._wp 
    326356      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 
    327361#if defined key_asminc 
    328362      ALLOCATE( ssh_iau      (jpi,jpj)     )   ;   ssh_iau      (:,:)   = 0._wp 
    … …  
    398432            ! to allow for differences in masks 
    399433            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 
    400570         ENDIF 
    401571         ! 
    … …  
    823993      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation 
    824994      ! 
    825       INTEGER  ::   it 
     995      INTEGER  ::   it, jl 
    826996      REAL(wp) ::   zincwgt   ! IAU weight for current time step 
    827997#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 
    8301004#endif 
    8311005      !!---------------------------------------------------------------------- 
    … …  
    8501024            ! 
    8511025#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 
    8701161            ! 
    8711162            ! seaice salinity balancing (to add) 
    … …  
    8951186            ! 
    8961187            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 
    9231194            ! 
    9241195#if defined key_cice && defined key_asminc 
    … …  
    9371208            ! 
    9381209         ENDIF 
    939  
    940 !#if defined defined key_si3 || defined key_cice 
    941 ! 
    942 !            IF (ln_seaicebal ) THEN        
    943 !             !! balancing salinity increments 
    944 !             !! simple case from limflx.F90 (doesn't include a mass flux) 
    945 !             !! assumption is that as ice concentration is reduced or increased 
    946 !             !! the snow and ice depths remain constant 
    947 !             !! note that snow is being created where ice concentration is being increased 
    948 !             !! - could be more sophisticated and 
    949 !             !! not do this (but would need to alter h_snow) 
    950 ! 
    951 !             usave(:,:,:)=sb(:,:,:)   ! use array as a temporary store 
    952 ! 
    953 !             DO jj = 1, jpj 
    954 !               DO ji = 1, jpi  
    955 !           ! calculate change in ice and snow mass per unit area 
    956 !           ! positive values imply adding salt to the ocean (results from ice formation) 
    957 !           ! fwf : ice formation and melting 
    958 ! 
    959 !                 zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt 
    960 ! 
    961 !           ! change salinity down to mixed layer depth 
    962 !                 mld=hmld_kara(ji,jj) 
    963 ! 
    964 !           ! prevent small mld 
    965 !           ! less than 10m can cause salinity instability  
    966 !                 IF (mld < 10) mld=10 
    967 ! 
    968 !           ! set to bottom of a level  
    969 !                 DO jk = jpk-1, 2, -1 
    970 !                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN  
    971 !                     mld=gdepw(ji,jj,jk+1) 
    972 !                     jkmax=jk 
    973 !                   ENDIF 
    974 !                 ENDDO 
    975 ! 
    976 !            ! avoid applying salinity balancing in shallow water or on land 
    977 !            !  
    978 ! 
    979 !            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m) 
    980 ! 
    981 !                 dsal_ocn=0.0_wp 
    982 !                 sal_thresh=5.0_wp        ! minimum salinity threshold for salinity balancing 
    983 ! 
    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 layer 
    988 !           ! but prevent salinity below a threshold value  
    989 ! 
    990 !                   DO jk = 1, jkmax               
    991 ! 
    992 !                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN  
    993 !                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn 
    994 !                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn 
    995 !                     ENDIF 
    996 ! 
    997 !                   ENDDO 
    998 ! 
    999 !      !            !  salt exchanges at the ice/ocean interface 
    1000 !      !            zpmess         = zfons / rdt_ice    ! rdt_ice is ice timestep 
    1001 !      ! 
    1002 !      !! Adjust fsalt. A +ve fsalt means adding salt to ocean 
    1003 !      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt   
    1004 !      !!                
    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 !ji 
    1009 !             ENDDO !jj! 
    1010 ! 
    1011 !            ENDIF !ln_seaicebal 
    1012 ! 
    1013 !#endif 
    10141210         ! 
    10151211      ENDIF 
    … …  
    10171213   END SUBROUTINE seaice_asm_inc 
    10181214    
     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 
    10191282   !!====================================================================== 
    10201283END MODULE asminc 
Note: See TracChangeset for help on using the changeset viewer.