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 15241 – NEMO

Changeset 15241


Ignore:
Timestamp:
2021-09-10T10:44:30+02:00 (3 years ago)
Author:
dcarneir
Message:

Cleaning the SIC IAU code and keeping only options used

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  
    9898#endif 
    9999   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 
    103101   LOGICAL, PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: incr_newice !: mask .TRUE.=DA positive ice increment to open water 
    104102#if defined key_cice && defined key_asminc 
     
    147145 
    148146#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 
    152147      REAL(wp) :: zopenwater_lim 
    153148      INTEGER  ::  jl 
     
    166161      ln_seaiceinc   = .FALSE. 
    167162      ln_temnofreeze = .FALSE. 
    168       na_iincr_split = 3 
    169       nh_iincr_min = 3 
    170163      zhi_damin = 0.5_wp 
    171164      zopenwater_lim = 1.0e-2_wp 
     
    199192         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix 
    200193         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 
    203194         WRITE(numout,*) '      Minimum ice thickness for new ice from da                zhi_damin = ', zhi_damin 
    204195         WRITE(numout,*) '      Limit for open water detection for ice da           zopenwater_lim = ', zopenwater_lim 
     
    441432            !convert single category increment to multi category 
    442433            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 
    467441               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 
    487467               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 
    570474         ENDIF 
    571475         ! 
     
    998902      REAL(wp), DIMENSION(jpi,jpj,jpl) :: da_i  ! change in ice concentration 
    999903      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 
    1001904      REAL(wp) :: rn_hinew_save 
    1002905      LOGICAL, SAVE :: initial_step=.TRUE. 
     
    1024927            ! 
    1025928#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 ) 
    1031930               z1_a_i(:,:,:) = 1.0_wp/a_i(:,:,:) 
    1032             elsewhere 
     931            ELSEWHERE 
    1033932               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 ) 
    1038948               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 
    1040951             
    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  
    1047952            ! compute change in ice concentration (new / old) 
    1048             where (incr_newice) 
     953            WHERE ( incr_newice ) 
    1049954               da_i(:,:,:) = 1.0_wp 
    1050             elsewhere 
     955            ELSEWHERE 
    1051956               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 
    11411017 
    11421018            ! make volume consistent with concentration and thickness 
     
    11501026            a_ip (:,:,:) = a_ip(:,:,:) * da_i(:,:,:) 
    11511027            v_ip (:,:,:) = v_ip(:,:,:) * da_i(:,:,:) 
    1152             do jl = 1, nlay_i 
     1028            DO jl = 1, nlay_i 
    11531029               e_i(:,:,jl,:) = e_i(:,:,jl,:) * da_i(:,:,:) 
    1154             end do 
    1155             do jl = 1, nlay_s 
     1030            END DO 
     1031            DO jl = 1, nlay_s 
    11561032               e_s(:,:,jl,:) = e_s(:,:,jl,:) * da_i(:,:,:) 
    1157             end do 
     1033            END DO 
    11581034 
    11591035            call ice_var_zapsmall()  
Note: See TracChangeset for help on using the changeset viewer.