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

Changeset 15476


Ignore:
Timestamp:
2021-11-08T10:31:34+01:00 (7 months ago)
Author:
dcarneir
Message:

Code review for IAU SIC with SI3

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.4_FOAM_IAU_SI3_SIC/src/OCE/ASM/asminc.F90

    r15412 r15476  
    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 
     37#if defined key_si3 && defined key_asminc 
    3838   USE phycst         ! physical constants 
    3939   USE ice1D          ! sea-ice: thermodynamics variables 
     
    9292   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment 
    9393   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc 
    94 #if defined key_si3 
     94   REAL(wp) :: zhi_damin                                            ! Ice thickness for new sea ice from DA increment 
     95#if defined key_si3 && defined key_asminc 
    9596   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   a_i_bkginc          ! Increment to the background sea ice conc categories 
    9697   LOGICAL, PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: incr_newice    ! Mask .TRUE.=DA positive ice increment to open water 
    97    REAL(wp) :: zhi_damin                                            ! Ice thickness for new sea ice from DA increment 
    9898#endif 
    9999#if defined key_cice && defined key_asminc 
     
    336336      ALLOCATE( ssh_bkginc   (jpi,jpj)     )   ;   ssh_bkginc   (:,:)   = 0._wp 
    337337      ALLOCATE( seaice_bkginc(jpi,jpj)     )   ;   seaice_bkginc(:,:)   = 0._wp 
    338 #if defined key_si3 
     338#if defined key_si3 && defined key_asminc 
    339339      ALLOCATE( a_i_bkginc   (jpi,jpj,jpl) )   ;   a_i_bkginc   (:,:,:) = 0._wp 
    340340      ALLOCATE( incr_newice(jpi,jpj,jpl) )     ;   incr_newice(:,:,:) = .FALSE. 
     
    414414            ! very small increments are also set to zero 
    415415            WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 .OR. & 
    416               &    ABS( seaice_bkginc(:,:) ) < epsi03 ) seaice_bkginc(:,:) = 0.0_wp 
     416              &    ABS( seaice_bkginc(:,:) ) < 1.0e-03_wp ) seaice_bkginc(:,:) = 0.0_wp 
    417417             
     418#if defined key_si3 && defined key_asminc 
    418419            IF (lwp) THEN 
    419420               WRITE(numout,*) 
    420                WRITE(numout,*) 'asm_inc_init : Converting seaice_bkginc to a_i_bkginc' 
     421               WRITE(numout,*) 'asm_inc_init : Converting single category increment to multi-category' 
    421422               WRITE(numout,*) '~~~~~~~~~~~~' 
    422423            END IF 
     
    424425            !single category increment for sea ice conc 
    425426            !convert single category increment to multi category 
    426             a_i_bkginc = 0.0_wp 
    427427            at_i = SUM( a_i(:,:,:), DIM=3 ) 
    428428 
    429429            ! ensure zhi_damin lies within 1st category 
    430             IF ( zhi_damin > hi_max(1) ) zhi_damin = hi_max(1) - epsi02 
     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 
    431438             
    432439            DO jj = 1, jpj 
    433440               DO ji = 1, jpi 
    434                   IF ( seaice_bkginc(ji,jj) > 0.0_wp) THEN 
     441                  IF ( seaice_bkginc(ji,jj) > 0.0_wp ) THEN 
    435442                     !positive ice concentration increments are always  
    436443                     !added to the thinnest category of ice 
     
    440447                     !available category until it reaches zero concentration 
    441448                     !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  
    442453                     zremaining_increment = seaice_bkginc(ji,jj) 
    443454                     DO jl = 1, jpl 
     
    451462            END DO 
    452463            ! find model points where DA new ice should be added to open water 
     464            ! in any ice category 
    453465            DO jl = 1,jpl 
    454466               WHERE ( at_i(:,:) < epsi02 .AND. seaice_bkginc(:,:) > 0.0_wp ) 
     
    457469            END DO 
    458470         ENDIF 
     471#endif 
    459472         ! 
    460473         CALL iom_close( inum ) 
     
    881894      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation 
    882895      ! 
    883       INTEGER  ::   it, jl 
     896      INTEGER  ::   it, jk 
    884897      REAL(wp) ::   zincwgt                            ! IAU weight for current time step 
    885 #if defined key_si3 
    886       LOGICAL, SAVE :: logical_newice = .TRUE.         ! Flag to add new ice from DA to open water 
     898#if defined key_si3 && defined key_asminc 
    887899      REAL(wp), DIMENSION(jpi,jpj,jpl) :: da_i         ! change in ice concentration 
    888900      REAL(wp), DIMENSION(jpi,jpj,jpl) :: dv_i         ! change in ice volume 
    889       REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i       ! inverse of old ice concentration 
    890       REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_v_i       ! inverse of old 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 
    891903      REAL(wp), DIMENSION(jpi,jpj)     :: zhi_damin_2D ! array with DA thickness for incr_newice 
    892904#endif 
     
    911923            ! Sea-ice : SI3 case 
    912924            ! 
    913 #if defined key_si3 
    914             ! ensure zhi_damin lies within 1st category 
    915             IF ( zhi_damin > hi_max(1) ) zhi_damin = hi_max(1) - epsi02 
    916  
     925#if defined key_si3 && defined key_asminc 
    917926            ! compute the inverse of key sea ice variables 
    918927            ! to be used later in the code  
     
    953962            ! just do it once since next IAU steps assume that new ice has  
    954963            ! already been added in 
    955             IF ( kt == nitiaustr_r .AND. logical_newice ) THEN 
     964            IF ( kt == nitiaustr_r ) THEN 
    956965 
    957966               ! assign zhi_damin to ice forming in open water 
     
    972981                  sv_i(:,:,:) = 0.0_wp 
    973982               END WHERE 
    974                DO jl = 1, nlay_i 
     983               DO jk = 1, nlay_i 
    975984                  WHERE ( incr_newice ) 
    976                      e_i(:,:,jl,:) = 0.0_wp 
     985                     e_i(:,:,jk,:) = 0.0_wp 
    977986                  END WHERE 
    978987               END DO 
    979                DO jl = 1, nlay_s 
     988               DO jk = 1, nlay_s 
    980989                  WHERE ( incr_newice ) 
    981                      e_s(:,:,jl,:) = 0.0_wp 
     990                     e_s(:,:,jk,:) = 0.0_wp 
    982991                  END WHERE 
    983992               END DO 
     
    985994               ! Initialisation of the salt content and ice enthalpy 
    986995               ! set flag of new ice to false after this 
    987                call init_new_ice_thd( zhi_damin_2D ) 
     996               CALL init_new_ice_thd( zhi_damin_2D ) 
    988997               incr_newice(:,:,:) = .FALSE. 
    989998            END IF 
     
    9911000            ! maintain equivalent values for key prognostic variables 
    9921001            v_s(:,:,:) = v_s(:,:,:) * da_i(:,:,:) 
    993             DO jl = 1, nlay_s 
    994                e_s(:,:,jl,:) = e_s(:,:,jl,:) * da_i(:,:,:) 
     1002            DO jk = 1, nlay_s 
     1003               e_s(:,:,jk,:) = e_s(:,:,jk,:) * da_i(:,:,:) 
    9951004            END DO 
    9961005            a_ip (:,:,:) = a_ip(:,:,:) * da_i(:,:,:) 
     
    9991008            ! ice volume dependent variables 
    10001009            sv_i (:,:,:) = sv_i(:,:,:) * dv_i(:,:,:) 
    1001             DO jl = 1, nlay_i 
    1002                e_i(:,:,jl,:) = e_i(:,:,jl,:) * dv_i(:,:,:) 
     1010            DO jk = 1, nlay_i 
     1011               e_i(:,:,jk,:) = e_i(:,:,jk,:) * dv_i(:,:,:) 
    10031012            END DO 
    10041013#endif 
     
    10111020            IF ( kt == nitiaufin_r ) THEN 
    10121021               DEALLOCATE( seaice_bkginc ) 
    1013 #if defined key_si3 
     1022#if defined key_si3 && defined key_asminc 
    10141023               DEALLOCATE( incr_newice ) 
    10151024               DEALLOCATE( a_i_bkginc ) 
     
    10661075      !!                forming at 1st category with thickness hi_new 
    10671076      !! 
    1068       !! ** Method  :   Apply SI3 thermodynamic equations to initialise 
    1069       !!                thermodynamic of new ice 
    1070       !! 
    1071       !! ** Action  :   update thermodynamic variables 
     1077      !! ** Method  :   Apply SI3 equations to initialise 
     1078      !!                thermodynamics of new ice 
     1079      !! 
     1080      !! ** Action  :   update sea ice thermodynamics 
    10721081      !!---------------------------------------------------------------------- 
    10731082      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    ::   hi_new  ! total thickness of new ice 
    10741083 
    1075       INTEGER :: jj, ji, jk, jl 
     1084      INTEGER :: jj, ji, jk 
    10761085      REAL(wp) ::   ztmelts   ! melting point 
    10771086 
     
    11081117         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)] 
    11091118            DO ji = 1, npti 
    1110                zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) ) 
     1119               zs_newice(ji) = MIN(  4.606_wp + 0.91_wp / zh_newice(ji) , rn_simax , 0.5_wp * sss_1d(ji) ) 
    11111120            END DO 
    11121121         CASE ( 3 )                    ! Sice = F(z) [multiyear ice] 
    1113             zs_newice(1:npti) =   2.3 
     1122            zs_newice(1:npti) =   2.3_wp 
    11141123         END SELECT 
    11151124 
     
    11231132         DO ji = 1, npti 
    11241133               ztmelts       = - rTmlt * zs_newice(ji)                  ! Melting point (C) 
    1125                e_i_1d(ji,:)  =   rhoi * (  rcpi  * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     & 
    1126                   &                      + rLfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   & 
     1134               e_i_1d(ji,:)  =   rhoi * (  rcpi  * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                        & 
     1135                  &                      + rLfus * ( 1.0_wp - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   & 
    11271136                  &                      - rcp   *         ztmelts ) 
    11281137         END DO 
Note: See TracChangeset for help on using the changeset viewer.