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

Changeset 15486


Ignore:
Timestamp:
2021-11-09T11:42:25+01:00 (7 months ago)
Author:
dcarneir
Message:

SIC assimilation interface (IAU) for SI3

Location:
NEMO/branches/UKMO/NEMO_4.0.4_FOAM_package/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.4_FOAM_package/src/ICE/ice.F90

    r14075 r15486  
    234234   REAL(wp), PUBLIC ::   rswitch          !: switch for the presence of ice (1) or not (0) 
    235235   REAL(wp), PUBLIC ::   rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft   !: conservation diagnostics 
     236   REAL(wp), PUBLIC, PARAMETER ::   epsi02 = 1.e-02_wp  !: small number  
     237   REAL(wp), PUBLIC, PARAMETER ::   epsi03 = 1.e-03_wp  !: small number  
    236238   REAL(wp), PUBLIC, PARAMETER ::   epsi06 = 1.e-06_wp  !: small number  
    237239   REAL(wp), PUBLIC, PARAMETER ::   epsi10 = 1.e-10_wp  !: small number  
  • NEMO/branches/UKMO/NEMO_4.0.4_FOAM_package/src/OCE/ASM/asminc.F90

    r14075 r15486  
    3535   USE sbc_oce         ! Surface boundary condition variables. 
    3636   USE diaobs   , ONLY : calc_date     ! Compute the calendar date on a given step 
    37 #if defined key_si3 
    38    USE ice      , ONLY : hm_i, at_i, at_i_b 
     37#if defined key_si3 && defined key_asminc 
     38   USE phycst         ! physical constants 
     39   USE ice1D          ! sea-ice: thermodynamics variables 
     40   USE icetab         ! sea-ice: 3D <==> 2D, 2D <==> 1D 
     41   USE ice 
    3942#endif 
    4043   ! 
     
    8992   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment 
    9093   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc 
     94   REAL(wp) :: zhi_damin                                            ! Ice thickness for new sea ice from DA increment 
     95#if defined key_si3 && defined key_asminc 
     96   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   a_i_bkginc          ! Increment to the background sea ice conc categories 
     97   LOGICAL, PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: incr_newice    ! Mask .TRUE.=DA positive ice increment to open water 
     98#endif 
    9199#if defined key_cice && defined key_asminc 
    92100   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ndaice_da             ! ice increment tendency into CICE 
     
    112120      !! ** Action  :  
    113121      !!---------------------------------------------------------------------- 
    114       INTEGER :: ji, jj, jk, jt  ! dummy loop indices 
    115       INTEGER :: imid, inum      ! local integers 
    116       INTEGER :: ios             ! Local integer output status for namelist read 
    117       INTEGER :: iiauper         ! Number of time steps in the IAU period 
    118       INTEGER :: icycper         ! Number of time steps in the cycle 
     122      INTEGER :: ji, jj, jk, jt, jl  ! dummy loop indices 
     123      INTEGER :: imid, inum          ! local integers 
     124      INTEGER :: ios                 ! Local integer output status for namelist read 
     125      INTEGER :: iiauper             ! Number of time steps in the IAU period 
     126      INTEGER :: icycper             ! Number of time steps in the cycle 
    119127      REAL(KIND=dp) :: ditend_date     ! Date YYYYMMDD.HHMMSS of final time step 
    120128      REAL(KIND=dp) :: ditbkg_date     ! Date YYYYMMDD.HHMMSS of background time step for Jb term 
     
    131139      ! 
    132140      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zhdiv   ! 2D workspace 
     141      REAL(wp) :: zremaining_increment 
     142 
    133143      !! 
    134144      NAMELIST/nam_asminc/ ln_bkgwri,                                      & 
    135145         &                 ln_trainc, ln_dyninc, ln_sshinc,                & 
    136          &                 ln_asmdin, ln_asmiau,                           & 
     146         &                 ln_seaiceinc, ln_asmdin, ln_asmiau,             & 
    137147         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   & 
    138          &                 ln_salfix, salfixmin, nn_divdmp 
     148         &                 ln_salfix, salfixmin, nn_divdmp, zhi_damin 
    139149      !!---------------------------------------------------------------------- 
    140150 
     
    173183         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix 
    174184         WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin 
     185         WRITE(numout,*) '      Minimum ice thickness for new ice from DA                zhi_damin = ', zhi_damin 
    175186      ENDIF 
    176187 
     
    325336      ALLOCATE( ssh_bkginc   (jpi,jpj)     )   ;   ssh_bkginc   (:,:)   = 0._wp 
    326337      ALLOCATE( seaice_bkginc(jpi,jpj)     )   ;   seaice_bkginc(:,:)   = 0._wp 
     338#if defined key_si3 && defined key_asminc 
     339      ALLOCATE( a_i_bkginc   (jpi,jpj,jpl) )   ;   a_i_bkginc   (:,:,:) = 0._wp 
     340      ALLOCATE( incr_newice(jpi,jpj,jpl) )     ;   incr_newice(:,:,:) = .FALSE. 
     341#endif 
    327342#if defined key_asminc 
    328343      ALLOCATE( ssh_iau      (jpi,jpj)     )   ;   ssh_iau      (:,:)   = 0._wp 
     
    397412            ! Set missing increments to 0.0 rather than 1e+20 
    398413            ! to allow for differences in masks 
    399             WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0 
    400          ENDIF 
     414            ! very small increments are also set to zero 
     415            WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 .OR. & 
     416              &    ABS( seaice_bkginc(:,:) ) < 1.0e-03_wp ) seaice_bkginc(:,:) = 0.0_wp 
     417             
     418#if defined key_si3 && defined key_asminc 
     419            IF (lwp) THEN 
     420               WRITE(numout,*) 
     421               WRITE(numout,*) 'asm_inc_init : Converting single category increment to multi-category' 
     422               WRITE(numout,*) '~~~~~~~~~~~~' 
     423            END IF 
     424 
     425            !single category increment for sea ice conc 
     426            !convert single category increment to multi category 
     427            at_i = SUM( a_i(:,:,:), DIM=3 ) 
     428 
     429            ! ensure zhi_damin lies within 1st category 
     430            IF ( zhi_damin > hi_max(1) ) THEN  
     431               IF (lwp) THEN 
     432                  WRITE(numout,*) 
     433                  WRITE(numout,*) 'minimum DA thickness > hmax(1): setting minimum DA thickness to hmax(1)' 
     434                  WRITE(numout,*) '~~~~~~~~~~~~' 
     435               END IF 
     436               zhi_damin = hi_max(1) - epsi02 
     437            END IF 
     438             
     439            DO jj = 1, jpj 
     440               DO ji = 1, jpi 
     441                  IF ( seaice_bkginc(ji,jj) > 0.0_wp ) THEN 
     442                     !positive ice concentration increments are always  
     443                     !added to the thinnest category of ice 
     444                     a_i_bkginc(ji,jj,1) = seaice_bkginc(ji,jj) 
     445                  ELSE 
     446                     !negative increments are first removed from the thinnest 
     447                     !available category until it reaches zero concentration 
     448                     !and then progressively removed from thicker categories 
     449 
     450                     !NOTE: might consider the remote possibility that zremaining_increment  
     451                     !left over is not zero, particulary if SI3 is being run  
     452                     !as a single category  
     453                     zremaining_increment = seaice_bkginc(ji,jj) 
     454                     DO jl = 1, jpl 
     455                        ! assign as much increment as possible to current category 
     456                        a_i_bkginc(ji,jj,jl) = -MIN( a_i(ji,jj,jl), -zremaining_increment ) 
     457                        ! update remaining amount of increment 
     458                        zremaining_increment = zremaining_increment - a_i_bkginc(ji,jj,jl) 
     459                     END DO 
     460                  END IF 
     461               END DO 
     462            END DO 
     463            ! find model points where DA new ice should be added to open water 
     464            ! in any ice category 
     465            DO jl = 1,jpl 
     466               WHERE ( at_i(:,:) < epsi02 .AND. seaice_bkginc(:,:) > 0.0_wp ) 
     467                  incr_newice(:,:,jl) = .TRUE. 
     468               END WHERE 
     469            END DO 
     470         ENDIF 
     471#endif 
    401472         ! 
    402473         CALL iom_close( inum ) 
     
    823894      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation 
    824895      ! 
    825       INTEGER  ::   it 
    826       REAL(wp) ::   zincwgt   ! IAU weight for current time step 
    827 #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 
     896      INTEGER  ::   it, jk 
     897      REAL(wp) ::   zincwgt                            ! IAU weight for current time step 
     898#if defined key_si3 && defined key_asminc 
     899      REAL(wp), DIMENSION(jpi,jpj,jpl) :: da_i         ! change in ice concentration 
     900      REAL(wp), DIMENSION(jpi,jpj,jpl) :: dv_i         ! change in ice volume 
     901      REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i       ! inverse of ice concentration before current IAU step 
     902      REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_v_i       ! inverse of ice volume before current IAU step 
     903      REAL(wp), DIMENSION(jpi,jpj)     :: zhi_damin_2D ! array with DA thickness for incr_newice 
    830904#endif 
    831905      !!---------------------------------------------------------------------- 
     
    849923            ! Sea-ice : SI3 case 
    850924            ! 
    851 #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     
     925#if defined key_si3 && defined key_asminc 
     926            ! compute the inverse of key sea ice variables 
     927            ! to be used later in the code  
     928            WHERE( a_i(:,:,:) > epsi10 ) 
     929               z1_a_i(:,:,:) = 1.0_wp / a_i(:,:,:) 
     930               z1_v_i(:,:,:) = 1.0_wp / v_i(:,:,:) 
    864931            ELSEWHERE 
    865                zhicifinc(:,:) = 0.0_wp 
     932               z1_a_i(:,:,:) = 0.0_wp 
     933               z1_v_i(:,:,:) = 0.0_wp 
    866934            END WHERE 
    867             ! 
    868             ! nudge ice depth 
    869             hm_i (:,:) = hm_i (:,:) + zhicifinc(:,:) 
    870             ! 
    871             ! seaice salinity balancing (to add) 
    872 #endif 
    873             ! 
     935 
     936            ! add positive concentration increments to regions where ice 
     937            ! is already present and bound them to 1 
     938            ! ice volume is added based on zhi_damin 
     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 ) 
     948               a_i(:,:,:) = MAX( a_i(:,:,:) + a_i_bkginc(:,:,:) * zincwgt, 0.0_wp ) 
     949               v_i(:,:,:) = a_i(:,:,:) * h_i(:,:,:) 
     950            END WHERE 
     951            
     952            ! compute changes in ice concentration and volume 
     953            WHERE ( incr_newice ) 
     954               da_i(:,:,:) = 1.0_wp 
     955               dv_i(:,:,:) = 1.0_wp 
     956            ELSEWHERE 
     957               da_i(:,:,:) = a_i(:,:,:) * z1_a_i(:,:,:) 
     958               dv_i(:,:,:) = v_i(:,:,:) * z1_v_i(:,:,:) 
     959            END WHERE 
     960 
     961            ! initialise thermodynamics of new ice being added to open water 
     962            ! just do it once since next IAU steps assume that new ice has  
     963            ! already been added in 
     964            IF ( kt == nitiaustr_r ) THEN 
     965 
     966               ! assign zhi_damin to ice forming in open water 
     967               WHERE ( ANY( incr_newice, DIM=3 ) ) 
     968                  zhi_damin_2D(:,:) = zhi_damin 
     969               ELSEWHERE 
     970                  zhi_damin_2D(:,:) = 0.0_wp 
     971               END WHERE 
     972 
     973               ! add ice concentration and volume 
     974               ! ensure the other prognostic variables are set to zero  
     975               WHERE ( incr_newice ) 
     976                  a_i(:,:,:) = MIN( 1.0_wp, a_i_bkginc(:,:,:) * zincwgt ) 
     977                  v_i(:,:,:) = MIN( 1.0_wp, a_i_bkginc(:,:,:) * zincwgt ) * zhi_damin 
     978                  v_s (:,:,:) = 0.0_wp 
     979                  a_ip(:,:,:) = 0.0_wp 
     980                  v_ip(:,:,:) = 0.0_wp 
     981                  sv_i(:,:,:) = 0.0_wp 
     982               END WHERE 
     983               DO jk = 1, nlay_i 
     984                  WHERE ( incr_newice ) 
     985                     e_i(:,:,jk,:) = 0.0_wp 
     986                  END WHERE 
     987               END DO 
     988               DO jk = 1, nlay_s 
     989                  WHERE ( incr_newice ) 
     990                     e_s(:,:,jk,:) = 0.0_wp 
     991                  END WHERE 
     992               END DO 
     993 
     994               ! Initialisation of the salt content and ice enthalpy 
     995               ! set flag of new ice to false after this 
     996               CALL init_new_ice_thd( zhi_damin_2D ) 
     997               incr_newice(:,:,:) = .FALSE. 
     998            END IF 
     999 
     1000            ! maintain equivalent values for key prognostic variables 
     1001            v_s(:,:,:) = v_s(:,:,:) * da_i(:,:,:) 
     1002            DO jk = 1, nlay_s 
     1003               e_s(:,:,jk,:) = e_s(:,:,jk,:) * da_i(:,:,:) 
     1004            END DO 
     1005            a_ip (:,:,:) = a_ip(:,:,:) * da_i(:,:,:) 
     1006            v_ip (:,:,:) = v_ip(:,:,:) * da_i(:,:,:) 
     1007  
     1008            ! ice volume dependent variables 
     1009            sv_i (:,:,:) = sv_i(:,:,:) * dv_i(:,:,:) 
     1010            DO jk = 1, nlay_i 
     1011               e_i(:,:,jk,:) = e_i(:,:,jk,:) * dv_i(:,:,:) 
     1012            END DO 
     1013#endif 
     1014 
    8741015#if defined key_cice && defined key_asminc 
    8751016            ! Sea-ice : CICE case. Pass ice increment tendency into CICE 
     
    8791020            IF ( kt == nitiaufin_r ) THEN 
    8801021               DEALLOCATE( seaice_bkginc ) 
     1022#if defined key_si3 && defined key_asminc 
     1023               DEALLOCATE( incr_newice ) 
     1024               DEALLOCATE( a_i_bkginc ) 
     1025#endif 
    8811026            ENDIF 
    8821027            ! 
     
    8951040            ! 
    8961041            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 
     1042 
     1043            IF(lwp) THEN 
     1044               WRITE(numout,*) 
     1045               WRITE(numout,*) 'seaice_asm_inc : sea ice direct initialization at time step = ', kt 
     1046               WRITE(numout,*) '~~~~~~~~~~~~' 
     1047            ENDIF 
    9231048            ! 
    9241049#if defined key_cice && defined key_asminc 
     
    9371062            ! 
    9381063         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 
    10141064         ! 
    10151065      ENDIF 
    10161066      ! 
    10171067   END SUBROUTINE seaice_asm_inc 
    1018     
     1068 
     1069 
     1070   SUBROUTINE init_new_ice_thd( hi_new ) 
     1071      !!---------------------------------------------------------------------- 
     1072      !!                  ***  ROUTINE init_new_ice_thd  *** 
     1073      !! 
     1074      !! ** Purpose :   Initialise thermodynamics of new ice  
     1075      !!                forming at 1st category with thickness hi_new 
     1076      !! 
     1077      !! ** Method  :   Apply SI3  equations to initialise 
     1078      !!                thermodynamics of new ice 
     1079      !! 
     1080      !! ** Action  :   update sea ice thermodynamics 
     1081      !!---------------------------------------------------------------------- 
     1082      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    ::   hi_new  ! total thickness of new ice 
     1083 
     1084      INTEGER :: jj, ji, jk 
     1085      REAL(wp) ::   ztmelts         ! melting point 
     1086      REAL(wp) ::   Sice_Fz=2.3_wp  ! Salinity of the ice = F(z) [multiyear ice] 
     1087 
     1088      REAL(wp), DIMENSION(jpij) ::   zh_newice   ! 1d version of hi_new 
     1089      REAL(wp), DIMENSION(jpij) ::   zs_newice   ! salinity of new ice 
     1090      !!---------------------------------------------------------------------- 
     1091 
     1092      ! Identify grid points where new ice forms 
     1093      npti = 0   ;   nptidx(:) = 0 
     1094      DO jj = 1, jpj 
     1095         DO ji = 1, jpi 
     1096            IF ( hi_new(ji,jj) > 0._wp ) THEN 
     1097               npti = npti + 1 
     1098               nptidx( npti ) = (jj - 1) * jpi + ji 
     1099            ENDIF 
     1100         END DO 
     1101      END DO 
     1102 
     1103      ! Move from 2-D to 1-D vectors 
     1104      IF ( npti > 0 ) THEN 
     1105         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 
     1106         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) ) 
     1107         CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , hi_new        ) 
     1108         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d    (1:npti) , sss_m         ) 
     1109         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d   (1:npti) , t_bo          ) 
     1110         DO jk = 1, nlay_i 
     1111            CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,1) ) 
     1112         END DO 
     1113 
     1114         ! --- Salinity of new ice --- !  
     1115         SELECT CASE ( nn_icesal ) 
     1116         CASE ( 1 )                    ! Sice = constant  
     1117            zs_newice(1:npti) = rn_icesal 
     1118         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)] 
     1119            DO ji = 1, npti 
     1120               zs_newice(ji) = MIN(  4.606_wp + 0.91_wp / zh_newice(ji) , rn_simax , 0.5_wp * sss_1d(ji) ) 
     1121            END DO 
     1122         CASE ( 3 )                    ! Sice = F(z) [multiyear ice] 
     1123            zs_newice(1:npti) = Sice_Fz 
     1124         END SELECT 
     1125 
     1126         ! --- Update ice salt content --- ! 
     1127         DO ji = 1, npti 
     1128            sv_i_2d(ji,1) = sv_i_2d(ji,1) + zs_newice(ji) * ( v_i_2d(ji,1) ) 
     1129         END DO 
     1130 
     1131         ! --- Heat content of new ice --- ! 
     1132         ! We assume that new ice is formed at the seawater freezing point 
     1133         DO ji = 1, npti 
     1134               ztmelts       = - rTmlt * zs_newice(ji)                  ! Melting point (C) 
     1135               e_i_1d(ji,:)  =   rhoi * (  rcpi  * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                        & 
     1136                  &                      + rLfus * ( 1.0_wp - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   & 
     1137                  &                      - rcp   *         ztmelts ) 
     1138         END DO 
     1139 
     1140         ! Change units for e_i 
     1141         DO jk = 1, nlay_i 
     1142            e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) * v_i_2d(1:npti,1) * r1_nlay_i  
     1143         END DO 
     1144 
     1145         ! Reforming full thermodynamic variables 
     1146         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 
     1147         DO jk = 1, nlay_i 
     1148               CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,1) ) 
     1149         END DO 
     1150      END IF 
     1151 
     1152   END SUBROUTINE init_new_ice_thd 
     1153 
     1154 
    10191155   !!====================================================================== 
    10201156END MODULE asminc 
Note: See TracChangeset for help on using the changeset viewer.