Ignore:
Timestamp:
2019-08-02T18:31:52+02:00 (15 months ago)
Author:
clem
Message:

rewrite iceistate.F90, add new variables (optional) for input fields, and prepare the code for a simplified ice restart

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/iceistate.F90

    r11317 r11397  
    2222   USE eosbn2         ! equation of state 
    2323   USE domvvl         ! Variable volume 
    24    USE ice            ! sea-ice variables 
    25    USE icevar         ! ice_var_salprof 
     24   USE ice            ! sea-ice: variables 
     25   USE ice1D          ! sea-ice: thermodynamics variables 
     26   USE icetab         ! sea-ice: 1D <==> 2D transformation 
     27   USE icevar         ! sea-ice: operations 
    2628   ! 
    2729   USE in_out_manager ! I/O manager 
     
    3638   PUBLIC   ice_istate        ! called by icestp.F90 
    3739   PUBLIC   ice_istate_init   ! called by icestp.F90 
    38  
    39    INTEGER , PARAMETER ::   jpfldi = 6           ! maximum number of files to read 
    40    INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness (m)    at T-point 
    41    INTEGER , PARAMETER ::   jp_hts = 2           ! index of snow thicknes (m)    at T-point 
    42    INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction (%) at T-point 
    43    INTEGER , PARAMETER ::   jp_tsu = 4           ! index of ice surface temp (K)    at T-point 
    44    INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temp at T-point 
    45    INTEGER , PARAMETER ::   jp_smi = 6           ! index of ice sali at T-point 
    46    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read) 
    4740   ! 
    4841   !                             !! ** namelist (namini) ** 
    49    LOGICAL  ::   ln_iceini        ! initialization or not 
    50    LOGICAL  ::   ln_iceini_file   ! Ice initialization state from 2D netcdf file 
     42   LOGICAL  ::   ln_iceini        ! Ice initialization or not 
     43   LOGICAL  ::   ln_iceini_file   ! Ice initialization from 2D netcdf file 
    5144   REAL(wp) ::   rn_thres_sst     ! threshold water temperature for initial sea ice 
    5245   REAL(wp) ::   rn_hts_ini_n     ! initial snow thickness in the north 
     
    5447   REAL(wp) ::   rn_hti_ini_n     ! initial ice thickness in the north 
    5548   REAL(wp) ::   rn_hti_ini_s     ! initial ice thickness in the south 
    56    REAL(wp) ::   rn_ati_ini_n     ! initial leads area in the north 
    57    REAL(wp) ::   rn_ati_ini_s     ! initial leads area in the south 
    58    REAL(wp) ::   rn_smi_ini_n     ! initial salinity  
    59    REAL(wp) ::   rn_smi_ini_s     ! initial salinity 
    60    REAL(wp) ::   rn_tmi_ini_n     ! initial temperature 
    61    REAL(wp) ::   rn_tmi_ini_s     ! initial temperature 
    62     
     49   REAL(wp) ::   rn_ati_ini_n     ! initial concentration in the north 
     50   REAL(wp) ::   rn_ati_ini_s     ! initial concentration in the south 
     51   REAL(wp) ::   rn_smi_ini_n     ! initial salinity in the north 
     52   REAL(wp) ::   rn_smi_ini_s     ! initial salinity in the south 
     53   REAL(wp) ::   rn_tmi_ini_n     ! initial ice temperature in the north 
     54   REAL(wp) ::   rn_tmi_ini_s     ! initial ice temperature in the south 
     55   REAL(wp) ::   rn_tms_ini_n     ! initial snw temperature in the north 
     56   REAL(wp) ::   rn_tms_ini_s     ! initial swn temperature in the south 
     57   REAL(wp) ::   rn_apd_ini_n     ! initial pond fraction in the north 
     58   REAL(wp) ::   rn_apd_ini_s     ! initial pond fraction in the south 
     59   REAL(wp) ::   rn_hpd_ini_n     ! initial pond depth in the north 
     60   REAL(wp) ::   rn_hpd_ini_s     ! initial pond depth in the south 
     61 
     62   !                              ! if ln_iceini_file = T 
     63   INTEGER , PARAMETER ::   jpfldi = 9           ! maximum number of files to read 
     64   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness    (m) 
     65   INTEGER , PARAMETER ::   jp_hts = 2           ! index of snw thickness    (m) 
     66   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction     (-) 
     67   INTEGER , PARAMETER ::   jp_tsu = 4           ! index of ice surface temp (K) 
     68   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temperature  (K) 
     69   INTEGER , PARAMETER ::   jp_smi = 6           ! index of ice salinity     (g/kg) 
     70   INTEGER , PARAMETER ::   jp_tms = 7           ! index of snw temperature  (K) 
     71   INTEGER , PARAMETER ::   jp_apd = 8           ! index of pnd fraction     (-) 
     72   INTEGER , PARAMETER ::   jp_hpd = 9           ! index of pnd depth        (m) 
     73   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read) 
     74   !    
    6375   !!---------------------------------------------------------------------- 
    6476   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
     
    6880CONTAINS 
    6981 
    70    SUBROUTINE ice_istate 
     82   SUBROUTINE ice_istate( kt ) 
    7183      !!------------------------------------------------------------------- 
    7284      !!                    ***  ROUTINE ice_istate  *** 
     
    8799      !! 
    88100      !! ** Notes   : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even 
    89       !!              where there is no ice (clem: I do not know why, is it mandatory?)  
     101      !!              where there is no ice 
    90102      !!-------------------------------------------------------------------- 
     103      INTEGER, INTENT(in) ::   kt   ! time step  
     104      !! 
    91105      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices 
    92       INTEGER  ::   i_hemis, i_fill, jl0   ! local integers 
    93       REAL(wp) ::   ztmelts, zdh 
    94       REAL(wp) ::   zarg, zV, zconv, zdv, zfac 
     106      REAL(wp) ::   ztmelts 
    95107      INTEGER , DIMENSION(4)           ::   itest 
    96108      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d 
    97109      REAL(wp), DIMENSION(jpi,jpj)     ::   zswitch    ! ice indicator 
    98       REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file 
    99       REAL(wp), DIMENSION(jpi,jpj)     ::   zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 
    100       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zh_i_ini , za_i_ini                        !data by cattegories to fill 
     110      REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, ztm_s_ini            !data from namelist or nc file 
     111      REAL(wp), DIMENSION(jpi,jpj)     ::   zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 
     112      REAL(wp), DIMENSION(jpi,jpj)     ::   zapnd_ini, zhpnd_ini                       !data from namelist or nc file 
     113      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !temporary arrays 
     114      !! 
     115      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d 
    101116      !-------------------------------------------------------------------- 
    102117 
     
    105120      IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
    106121 
    107       !-------------------------------------------------------------------- 
    108       ! 1) Set surface and bottom temperatures to initial values 
    109       !-------------------------------------------------------------------- 
    110       ! 
    111       ! init surface temperature 
     122      !--------------------------- 
     123      ! 1) 1st init. of the fields 
     124      !--------------------------- 
     125      ! 
     126      ! basal temperature (considered at freezing point)   [Kelvin] 
     127      CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 
     128      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)  
     129      ! 
     130      ! surface temperature and conductivity 
    112131      DO jl = 1, jpl 
    113132         t_su   (:,:,jl) = rt0 * tmask(:,:,1)  ! temp at the surface 
     
    115134      END DO 
    116135      ! 
    117       ! init basal temperature (considered at freezing point)   [Kelvin] 
    118       CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 
    119       t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)  
    120  
     136      ! ice and snw temperatures 
     137      DO jl = 1, jpl 
     138         DO jk = 1, nlay_i 
     139            t_i(:,:,jk,jl) = rt0 * tmask(:,:,1) 
     140         END DO 
     141         DO jk = 1, nlay_s 
     142            t_s(:,:,jk,jl) = rt0 * tmask(:,:,1) 
     143         END DO 
     144      END DO 
     145      ! 
     146      ! specific temperatures for coupled runs 
     147      tn_ice (:,:,:) = t_i (:,:,1,:) 
     148      t1_ice (:,:,:) = t_i (:,:,1,:) 
     149 
     150      ! heat contents 
     151      e_i (:,:,:,:) = 0._wp 
     152      e_s (:,:,:,:) = 0._wp 
     153       
     154      ! general fields 
     155      a_i (:,:,:) = 0._wp 
     156      v_i (:,:,:) = 0._wp 
     157      v_s (:,:,:) = 0._wp 
     158      sv_i(:,:,:) = 0._wp 
     159      oa_i(:,:,:) = 0._wp 
     160      ! 
     161      h_i (:,:,:) = 0._wp 
     162      h_s (:,:,:) = 0._wp 
     163      s_i (:,:,:) = 0._wp 
     164      o_i (:,:,:) = 0._wp 
     165      ! 
     166      ! melt ponds 
     167      a_ip     (:,:,:) = 0._wp 
     168      v_ip     (:,:,:) = 0._wp 
     169      a_ip_frac(:,:,:) = 0._wp 
     170      h_ip     (:,:,:) = 0._wp 
     171      ! 
     172      ! ice velocities 
     173      u_ice (:,:) = 0._wp 
     174      v_ice (:,:) = 0._wp 
     175      ! 
     176      !------------------------------------------------------------------------ 
     177      ! 2) overwrite some of the fields with namelist parameters or netcdf file 
     178      !------------------------------------------------------------------------ 
    121179      IF( ln_iceini ) THEN 
    122          !----------------------------------------------------------- 
    123          ! 2) Compute or read sea ice variables ===> single category 
    124          !----------------------------------------------------------- 
    125          ! 
    126180         !                             !---------------! 
    127181         IF( ln_iceini_file )THEN      ! Read a file   ! 
    128182            !                          !---------------! 
    129             ! 
    130             zht_i_ini(:,:)  = si(jp_hti)%fnow(:,:,1) 
    131             zht_s_ini(:,:)  = si(jp_hts)%fnow(:,:,1) 
    132             zat_i_ini(:,:)  = si(jp_ati)%fnow(:,:,1) 
    133             zts_u_ini(:,:)  = si(jp_tsu)%fnow(:,:,1) 
    134             ztm_i_ini(:,:)  = si(jp_tmi)%fnow(:,:,1) 
    135             zsm_i_ini(:,:)  = si(jp_smi)%fnow(:,:,1) 
     183            CALL fld_read( kt, 1, si )                ! input fields provided at the current time-step 
     184            ! 
     185            zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) 
     186            zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) 
     187            zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) 
     188            zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) 
     189            ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) 
     190            zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) 
     191            ! snow temperature 
     192            IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' )   si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 
     193            ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) 
     194            ! ponds 
     195            IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' )   si(jp_apd)%fnow(:,:,1) = 0._wp  
     196            IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' )   si(jp_hpd)%fnow(:,:,1) = 0._wp  
     197            zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) 
     198            zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) 
    136199            ! 
    137200            WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1)  
    138201            ELSEWHERE                       ; zswitch(:,:) = 0._wp 
    139202            END WHERE 
    140             zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:) 
    141             ! 
    142203            !                          !---------------! 
    143204         ELSE                          ! Read namelist ! 
    144205            !                          !---------------! 
    145             ! no ice if sst <= t-freez + ttest 
     206            ! no ice if (sst - Tfreez) >= thresold 
    146207            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp  
    147208            ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1) 
     
    153214               zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:) 
    154215               zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:) 
    155                zts_u_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 
     216               zt_su_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 
    156217               zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:) 
    157218               ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 
     219               ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:) 
     220               zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) 
     221               zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:) 
    158222            ELSEWHERE 
    159223               zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:) 
    160224               zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:) 
    161225               zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:) 
    162                zts_u_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 
     226               zt_su_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 
    163227               zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:) 
    164228               ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 
     229               ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:) 
     230               zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) 
     231               zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:) 
    165232            END WHERE 
    166             zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:) 
    167233            ! 
    168234         ENDIF 
     235         !-------------! 
     236         ! fill fields ! 
     237         !-------------! 
     238         ! select ice covered grid points 
     239         npti = 0 ; nptidx(:) = 0 
     240         DO jj = 1, jpj 
     241            DO ji = 1, jpi 
     242               IF ( zht_i_ini(ji,jj) > 0._wp ) THEN 
     243                  npti         = npti  + 1 
     244                  nptidx(npti) = (jj - 1) * jpi + ji 
     245               ENDIF 
     246            END DO 
     247         END DO 
     248 
     249         ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj) 
     250         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti)  , zht_i_ini ) 
     251         CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti)  , zht_s_ini ) 
     252         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)  , zat_i_ini ) 
     253         CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini ) 
     254         CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini ) 
     255         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti)  , zt_su_ini ) 
     256         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti)  , zsm_i_ini ) 
     257 
     258         ! allocate temporary arrays 
     259         ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), & 
     260            &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl) ) 
    169261          
    170          !------------------------------------------------------------------ 
    171          ! 3) Distribute ice concentration and thickness into the categories 
    172          !------------------------------------------------------------------ 
    173          ! a gaussian distribution for ice concentration is used 
    174          ! then we check whether the distribution fullfills 
    175          ! volume and area conservation, positivity and ice categories bounds 
    176  
    177          IF( jpl == 1 ) THEN 
    178             ! 
    179             zh_i_ini(:,:,1) = zht_i_ini(:,:) 
    180             za_i_ini(:,:,1) = zat_i_ini(:,:)             
    181             ! 
    182          ELSE 
    183             zh_i_ini(:,:,:) = 0._wp  
    184             za_i_ini(:,:,:) = 0._wp 
    185             ! 
     262         ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl) 
     263         CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                   zhi_2d, zhs_2d, zai_2d , & 
     264            &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti),   zti_2d, zts_2d, ztsu_2d, zsi_2d ) 
     265 
     266         ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl) 
     267         DO jl = 1, jpl 
     268            zti_3d(:,:,jl) = rt0 * tmask(:,:,1) 
     269            zts_3d(:,:,jl) = rt0 * tmask(:,:,1) 
     270         END DO 
     271         CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d  , h_i  ) 
     272         CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d  , h_s  ) 
     273         CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d  , a_i  ) 
     274         CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d  , zti_3d ) 
     275         CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d  , zts_3d ) 
     276         CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su ) 
     277         CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d  , s_i  ) 
     278 
     279         ! deallocate temporary arrays 
     280         DEALLOCATE( zhi_2d, zhs_2d, zai_2d , & 
     281            &        zti_2d, zts_2d, ztsu_2d, zsi_2d ) 
     282 
     283         ! Melt ponds: distribute uniformely over the categories 
     284         IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN 
     285            DO jl = 1, jpl 
     286               a_ip_frac(:,:,jl) = zapnd_ini(:,:) 
     287               h_ip     (:,:,jl) = zhpnd_ini(:,:) 
     288               a_ip     (:,:,jl) = a_ip_frac(:,:,jl) * a_i (:,:,jl)  
     289               v_ip     (:,:,jl) = h_ip     (:,:,jl) * a_ip(:,:,jl) 
     290            END DO 
     291         ENDIF 
     292           
     293         ! calculate extensive and intensive variables 
     294         CALL ice_var_salprof ! for sz_i 
     295         DO jl = 1, jpl 
    186296            DO jj = 1, jpj 
    187297               DO ji = 1, jpi 
    188                   ! 
    189                   IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN 
    190  
    191                      ! find which category (jl0) the input ice thickness falls into 
    192                      jl0 = jpl 
    193                      DO jl = 1, jpl 
    194                         IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN 
    195                            jl0 = jl 
    196                            CYCLE 
    197                         ENDIF 
    198                      END DO 
    199                      ! 
    200                      itest(:) = 0 
    201                      i_fill   = jpl + 1                                            !------------------------------------ 
    202                      DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories 
    203                         !                                                          !------------------------------------ 
    204                         i_fill = i_fill - 1 
    205                         ! 
    206                         zh_i_ini(ji,jj,:) = 0._wp  
    207                         za_i_ini(ji,jj,:) = 0._wp 
    208                         itest(:) = 0 
    209                         ! 
    210                         IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1 
    211                            zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj) 
    212                            za_i_ini(ji,jj,1) = zat_i_ini(ji,jj) 
    213                         ELSE                         !-- case ice is thicker: fill categories >1 
    214                            ! thickness 
    215                            DO jl = 1, i_fill-1 
    216                               zh_i_ini(ji,jj,jl) = hi_mean(jl) 
    217                            END DO 
    218                            ! 
    219                            ! concentration 
    220                            za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 
    221                            DO jl = 1, i_fill - 1 
    222                               IF( jl /= jl0 )THEN 
    223                                  zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) ) 
    224                                  za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 
    225                               ENDIF 
    226                            END DO 
    227  
    228                            ! last category 
    229                            za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) ) 
    230                            zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) ) 
    231                            zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 )  
    232  
    233                            ! correction if concentration of upper cat is greater than lower cat 
    234                            !   (it should be a gaussian around jl0 but sometimes it is not) 
    235                            IF ( jl0 /= jpl ) THEN 
    236                               DO jl = jpl, jl0+1, -1 
    237                                  IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN 
    238                                     zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl) 
    239                                     zh_i_ini(ji,jj,jl    ) = 0._wp 
    240                                     za_i_ini(ji,jj,jl    ) = 0._wp 
    241                                     za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  & 
    242                                        &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 ) 
    243                                  END IF 
    244                               ENDDO 
    245                            ENDIF 
    246                            ! 
    247                         ENDIF 
    248                         ! 
    249                         ! Compatibility tests 
    250                         zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )           ! Test 1: area conservation 
    251                         IF ( zconv < epsi06 ) itest(1) = 1 
    252                         ! 
    253                         zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &         ! Test 2: volume conservation 
    254                            &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) ) 
    255                         IF ( zconv < epsi06 ) itest(2) = 1 
    256                         ! 
    257                         IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1           ! Test 3: thickness of the last category is in-bounds ? 
    258                         ! 
    259                         itest(4) = 1 
    260                         DO jl = 1, i_fill 
    261                            IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0                        ! Test 4: positivity of ice concentrations 
    262                         END DO 
    263                         !                                                          !---------------------------- 
    264                      END DO                                                        ! end iteration on categories 
    265                      !                                                             !---------------------------- 
    266                      IF( lwp .AND. SUM(itest) /= 4 ) THEN  
    267                         WRITE(numout,*) 
    268                         WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! ' 
    269                         WRITE(numout,*) ' !!!! Something is wrong in the SI3 initialization procedure ' 
    270                         WRITE(numout,*) 
    271                         WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:) 
    272                         WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 
    273                         WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 
    274                      ENDIF 
    275                      ! 
    276                   ENDIF 
    277                   ! 
    278                END DO 
    279             END DO 
    280          ENDIF 
    281           
    282          !--------------------------------------------------------------------- 
    283          ! 4) Fill in sea ice arrays 
    284          !--------------------------------------------------------------------- 
    285          ! 
    286          ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature 
    287          DO jl = 1, jpl ! loop over categories 
    288             DO jj = 1, jpj 
    289                DO ji = 1, jpi 
    290                   a_i(ji,jj,jl)  = zswitch(ji,jj) * za_i_ini(ji,jj,jl)                       ! concentration 
    291                   h_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(ji,jj,jl)                       ! ice thickness 
    292                   s_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(ji,jj)                         ! salinity 
    293                   o_i(ji,jj,jl)  = 0._wp                                                     ! age (0 day) 
    294                   t_su(ji,jj,jl) = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp 
    295                   ! 
    296                   IF( zht_i_ini(ji,jj) > 0._wp )THEN 
    297                     h_s(ji,jj,jl)= h_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) )  ! snow depth 
    298                   ELSE 
    299                     h_s(ji,jj,jl)= 0._wp 
    300                   ENDIF 
    301                   ! 
    302                   ! This case below should not be used if (h_s/h_i) is ok in namelist 
    303                   ! In case snow load is in excess that would lead to transformation from snow to ice 
    304                   ! Then, transfer the snow excess into the ice (different from icethd_dh) 
    305                   zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 )  
    306                   ! recompute h_i, h_s avoiding out of bounds values 
    307                   h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh ) 
    308                   h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi * r1_rhos ) 
    309                   ! 
    310                   ! ice volume, salt content, age content 
    311298                  v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume 
    312299                  v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume 
    313                   sv_i(ji,jj,jl) = MIN( s_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 
    314                   oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content 
     300                  sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl) 
    315301               END DO 
    316302            END DO 
    317303         END DO 
    318304         ! 
    319          IF( nn_icesal /= 2 )  THEN         ! for constant salinity in time 
    320             CALL ice_var_salprof 
    321             sv_i = s_i * v_i 
    322          ENDIF 
    323          !   
    324          ! Snow temperature and heat content 
    325          DO jk = 1, nlay_s 
    326             DO jl = 1, jpl ! loop over categories 
     305         DO jl = 1, jpl 
     306            DO jk = 1, nlay_s 
    327307               DO jj = 1, jpj 
    328308                  DO ji = 1, jpi 
    329                      t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 
    330                      ! Snow energy of melting 
    331                      e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) 
    332                      ! 
    333                      ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 
    334                      e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 
     309                     t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl) 
     310                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * & 
     311                        &               rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) 
    335312                  END DO 
    336313               END DO 
     
    338315         END DO 
    339316         ! 
    340          ! Ice salinity, temperature and heat content 
    341          DO jk = 1, nlay_i 
    342             DO jl = 1, jpl ! loop over categories 
     317         DO jl = 1, jpl 
     318            DO jk = 1, nlay_i 
    343319               DO jj = 1, jpj 
    344320                  DO ji = 1, jpi 
    345                      t_i (ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0  
    346                      sz_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin 
     321                     t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl)  
    347322                     ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 
    348                      ! 
    349                      ! heat content per unit volume 
    350                      e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoi * (   rcpi    * ( ztmelts - t_i(ji,jj,jk,jl) )           & 
    351                         &             + rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0) , -epsi20 )  )   & 
    352                         &             - rcp  * ( ztmelts - rt0 ) ) 
    353                      ! 
    354                      ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
    355                      e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 
     323                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * & 
     324                        &               rhoi * (  rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) + & 
     325                        &                         rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) & 
     326                        &                       - rcp   * ( ztmelts - rt0 ) ) 
    356327                  END DO 
    357328               END DO 
    358329            END DO 
    359330         END DO 
    360          ! 
    361          tn_ice (:,:,:) = t_su (:,:,:) 
    362          t1_ice (:,:,:) = t_i (:,:,1,:)   ! initialisation of 1st layer temp for coupled simu 
    363  
    364          ! Melt pond volume and fraction 
    365          IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN   ;   zfac = 1._wp 
    366          ELSE                                     ;   zfac = 0._wp 
    367          ENDIF  
    368          DO jl = 1, jpl 
    369             a_ip_frac(:,:,jl) = rn_apnd * zswitch(:,:) * zfac 
    370             h_ip     (:,:,jl) = rn_hpnd * zswitch(:,:) * zfac 
    371          END DO 
    372          a_ip(:,:,:) = a_ip_frac(:,:,:) * a_i (:,:,:)  
    373          v_ip(:,:,:) = h_ip     (:,:,:) * a_ip(:,:,:) 
    374          ! 
    375       ELSE ! if ln_iceini=false 
    376          a_i  (:,:,:) = 0._wp 
    377          v_i  (:,:,:) = 0._wp 
    378          v_s  (:,:,:) = 0._wp 
    379          sv_i (:,:,:) = 0._wp 
    380          oa_i (:,:,:) = 0._wp 
    381          h_i  (:,:,:) = 0._wp 
    382          h_s  (:,:,:) = 0._wp 
    383          s_i  (:,:,:) = 0._wp 
    384          o_i  (:,:,:) = 0._wp 
    385          ! 
    386          e_i(:,:,:,:) = 0._wp 
    387          e_s(:,:,:,:) = 0._wp 
    388          ! 
    389          DO jl = 1, jpl 
    390             DO jk = 1, nlay_i 
    391                t_i(:,:,jk,jl) = rt0 * tmask(:,:,1) 
    392             END DO 
    393             DO jk = 1, nlay_s 
    394                t_s(:,:,jk,jl) = rt0 * tmask(:,:,1) 
    395             END DO 
    396          END DO 
    397  
    398          tn_ice (:,:,:) = t_i (:,:,1,:) 
    399          t1_ice (:,:,:) = t_i (:,:,1,:)   ! initialisation of 1st layer temp for coupled simu 
    400           
    401          a_ip(:,:,:)      = 0._wp 
    402          v_ip(:,:,:)      = 0._wp 
    403          a_ip_frac(:,:,:) = 0._wp 
    404          h_ip     (:,:,:) = 0._wp 
     331 
     332         ! specific temperatures for coupled runs 
     333         tn_ice(:,:,:) = t_su(:,:,:) 
     334         t1_ice(:,:,:) = t_i (:,:,1,:) 
    405335         ! 
    406336      ENDIF ! ln_iceini 
    407337      ! 
    408       at_i (:,:) = 0.0_wp 
    409       DO jl = 1, jpl 
    410          at_i (:,:) = at_i (:,:) + a_i (:,:,jl) 
    411       END DO 
    412       ! 
    413       ! --- set ice velocities --- ! 
    414       u_ice (:,:) = 0._wp 
    415       v_ice (:,:) = 0._wp 
    416       ! fields needed for ice_dyn_adv_umx 
    417       l_split_advumx(1) = .FALSE. 
     338      at_i(:,:) = SUM( a_i, dim=3 ) 
    418339      ! 
    419340      !---------------------------------------------- 
    420       ! 5) Snow-ice mass (case ice is fully embedded) 
     341      ! 3) Snow-ice mass (case ice is fully embedded) 
    421342      !---------------------------------------------- 
    422343      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3  )   ! snow+ice mass 
     
    470391       
    471392      !------------------------------------ 
    472       ! 6) store fields at before time-step 
     393      ! 4) store fields at before time-step 
    473394      !------------------------------------ 
    474395      ! it is only necessary for the 1st interpolation by Agrif 
     
    508429      ! 
    509430      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files 
    510       TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi 
     431      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd 
    511432      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read 
    512433      ! 
    513       NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s,  & 
    514          &             rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, & 
    515          &             rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s,                             & 
    516          &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir 
     434      NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, & 
     435         &             rn_hts_ini_n, rn_hts_ini_s, rn_hti_ini_n, rn_hti_ini_s, & 
     436         &             rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, & 
     437         &             rn_tmi_ini_n, rn_tmi_ini_s, rn_tms_ini_n, rn_tms_ini_s, & 
     438         &             rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, & 
     439         &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir 
    517440      !!----------------------------------------------------------------------------- 
    518441      ! 
     
    528451      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu 
    529452      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi 
     453      slf_i(jp_tms) = sn_tms  ; 
     454      slf_i(jp_apd) = sn_apd  ;  slf_i(jp_hpd) = sn_hpd 
    530455      ! 
    531456      IF(lwp) THEN                          ! control print 
     
    534459         WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    535460         WRITE(numout,*) '   Namelist namini:' 
    536          WRITE(numout,*) '      initialization with ice (T) or not (F)                 ln_iceini       = ', ln_iceini 
    537          WRITE(numout,*) '      ice initialization from a netcdf file                  ln_iceini_file  = ', ln_iceini_file 
    538          WRITE(numout,*) '      max delta ocean temp. above Tfreeze with initial ice   rn_thres_sst    = ', rn_thres_sst 
    539          WRITE(numout,*) '      initial snow thickness in the north                    rn_hts_ini_n    = ', rn_hts_ini_n 
    540          WRITE(numout,*) '      initial snow thickness in the south                    rn_hts_ini_s    = ', rn_hts_ini_s  
    541          WRITE(numout,*) '      initial ice thickness  in the north                    rn_hti_ini_n    = ', rn_hti_ini_n 
    542          WRITE(numout,*) '      initial ice thickness  in the south                    rn_hti_ini_s    = ', rn_hti_ini_s 
    543          WRITE(numout,*) '      initial ice concentr.  in the north                    rn_ati_ini_n    = ', rn_ati_ini_n 
    544          WRITE(numout,*) '      initial ice concentr.  in the north                    rn_ati_ini_s    = ', rn_ati_ini_s 
    545          WRITE(numout,*) '      initial  ice salinity  in the north                    rn_smi_ini_n    = ', rn_smi_ini_n 
    546          WRITE(numout,*) '      initial  ice salinity  in the south                    rn_smi_ini_s    = ', rn_smi_ini_s 
    547          WRITE(numout,*) '      initial  ice/snw temp  in the north                    rn_tmi_ini_n    = ', rn_tmi_ini_n 
    548          WRITE(numout,*) '      initial  ice/snw temp  in the south                    rn_tmi_ini_s    = ', rn_tmi_ini_s 
     461         WRITE(numout,*) '      ice initialization (T) or not (F)                     ln_iceini       = ', ln_iceini 
     462         WRITE(numout,*) '      ice initialization from a netcdf file                 ln_iceini_file  = ', ln_iceini_file 
     463         WRITE(numout,*) '      max delta ocean temp. above Tfreeze with initial ice  rn_thres_sst    = ', rn_thres_sst 
     464         WRITE(numout,*) '      initial snw thickness in the north                    rn_hts_ini_n    = ', rn_hts_ini_n 
     465         WRITE(numout,*) '      initial snw thickness in the south                    rn_hts_ini_s    = ', rn_hts_ini_s  
     466         WRITE(numout,*) '      initial ice thickness in the north                    rn_hti_ini_n    = ', rn_hti_ini_n 
     467         WRITE(numout,*) '      initial ice thickness in the south                    rn_hti_ini_s    = ', rn_hti_ini_s 
     468         WRITE(numout,*) '      initial ice concentr  in the north                    rn_ati_ini_n    = ', rn_ati_ini_n 
     469         WRITE(numout,*) '      initial ice concentr  in the north                    rn_ati_ini_s    = ', rn_ati_ini_s 
     470         WRITE(numout,*) '      initial ice salinity  in the north                    rn_smi_ini_n    = ', rn_smi_ini_n 
     471         WRITE(numout,*) '      initial ice salinity  in the south                    rn_smi_ini_s    = ', rn_smi_ini_s 
     472         WRITE(numout,*) '      initial ice temperat  in the north                    rn_tmi_ini_n    = ', rn_tmi_ini_n 
     473         WRITE(numout,*) '      initial ice temperat  in the south                    rn_tmi_ini_s    = ', rn_tmi_ini_s 
     474         WRITE(numout,*) '      initial snw temperat  in the north                    rn_tms_ini_n    = ', rn_tms_ini_n 
     475         WRITE(numout,*) '      initial snw temperat  in the south                    rn_tms_ini_s    = ', rn_tms_ini_s 
     476         WRITE(numout,*) '      initial pnd fraction  in the north                    rn_apd_ini_n    = ', rn_apd_ini_n 
     477         WRITE(numout,*) '      initial pnd fraction  in the south                    rn_apd_ini_s    = ', rn_apd_ini_s 
     478         WRITE(numout,*) '      initial pnd depth     in the north                    rn_hpd_ini_n    = ', rn_hpd_ini_n 
     479         WRITE(numout,*) '      initial pnd depth     in the south                    rn_hpd_ini_s    = ', rn_hpd_ini_s 
    549480      ENDIF 
    550481      ! 
     
    554485         ALLOCATE( si(jpfldi), STAT=ierror ) 
    555486         IF( ierror > 0 ) THEN 
    556             CALL ctl_stop( 'Ice_ini in iceistate: unable to allocate si structure' )   ;   RETURN 
     487            CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' )   ;   RETURN 
    557488         ENDIF 
    558489         ! 
    559490         DO ifpr = 1, jpfldi 
    560491            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) ) 
    561             ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) ) 
     492            IF( slf_i(ifpr)%ln_tint )  ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) ) 
    562493         END DO 
    563494         ! 
    564495         ! fill si with slf_i and control print 
    565          CALL fld_fill( si, slf_i, cn_dir, 'ice_istate', 'ice istate ini', 'numnam_ice' ) 
    566          ! 
    567          CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step 
     496         CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' ) 
    568497         ! 
    569498      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.