Changeset 5976


Ignore:
Timestamp:
2015-12-02T12:46:21+01:00 (5 years ago)
Author:
cbricaud
Message:

code and namelist changes for ticket #1455: enable LIM3 initialisation with a netcdf file

Location:
branches/2015/dev_5974_MERCATOR10_LIMINI/NEMOGCM
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_5974_MERCATOR10_LIMINI/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r5429 r5976  
    3232!------------------------------------------------------------------------------ 
    3333   ln_iceini      = .true.         !  activate ice initialization (T) or not (F) 
     34   ln_iceini_file = .true.         !  ice initialization from a netcdf file (T) or constant (F) 
     35! 
     36!         !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
     37!         !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      ! 
     38   sn_hti = 'Ice_initialization'   , -12    ,'hti'      ,  .false.     , .true. , 'yearly'  , ''       , ''       , '' 
     39   sn_hts = 'Ice_initialization'   , -12    ,'hts'      ,  .false.     , .true. , 'yearly'  , ''       , ''       , '' 
     40   sn_ati = 'Ice_initialization'   , -12    ,'ati'      ,  .false.     , .true. , 'yearly'  , ''       , ''       , '' 
     41   sn_tsu = 'Ice_initialization'   , -12    ,'tsu'      ,  .false.     , .true. , 'yearly'  , ''       , ''       , '' 
     42   sn_tmi = 'Ice_initialization'   , -12    ,'tmi'      ,  .false.     , .true. , 'yearly'  , ''       , ''       , '' 
     43   sn_smi = 'Ice_initialization'   , -12    ,'smi'      ,  .false.     , .true. , 'yearly'  , ''       , ''       , '' 
     44   cn_dir='./' 
     45   ! 
    3446   rn_thres_sst   =  2.0           !  maximum water temperature with initial ice (degC) 
    3547   rn_hts_ini_n   =  0.3           !  initial real snow thickness (m), North 
  • branches/2015/dev_5974_MERCATOR10_LIMINI/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r5541 r5976  
    2828   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2929   USE wrk_nemo         ! work arrays 
     30   USE fldread          ! read input fields 
     31   USE iom 
    3032 
    3133   IMPLICIT NONE 
     
    4749   REAL(wp) ::   rn_tmi_ini_s   ! initial temperature 
    4850 
    49    LOGICAL  ::  ln_iceini    ! initialization or not 
     51   INTEGER , PARAMETER ::   jpfldi = 6           ! maximum number of files to read 
     52   INTEGER , PARAMETER ::   jp_hti = 1           ! index of thick (m)    at T-point 
     53   INTEGER , PARAMETER ::   jp_hts = 2           ! index of thick (m)    at T-point 
     54   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction (%) at T-point 
     55   INTEGER , PARAMETER ::   jp_tsu = 4           ! index of ice surface temp (K)    at T-point 
     56   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temp at T-point 
     57   INTEGER , PARAMETER ::   jp_smi = 6           ! index of ice sali at T-point 
     58   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si    ! structure of input fields (file informations, fields read) 
     59 
     60   LOGICAL  ::  ln_iceini        ! initialization or not 
     61   LOGICAL  ::  ln_iceini_file   ! Ice initialization state from 2D netcdf file 
    5062   !!---------------------------------------------------------------------- 
    5163   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008) 
     
    88100      REAL(wp)   :: ztmelts, zdh 
    89101      INTEGER    :: i_hemis, i_fill, jl0   
    90       REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 
    91       REAL(wp), POINTER, DIMENSION(:)     :: zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini 
    92       REAL(wp), POINTER, DIMENSION(:,:)   :: zh_i_ini, za_i_ini, zv_i_ini 
    93       REAL(wp), POINTER, DIMENSION(:,:)   :: zswitch    ! ice indicator 
    94       INTEGER,  POINTER, DIMENSION(:,:)   :: zhemis   ! hemispheric index 
    95       !-------------------------------------------------------------------- 
    96  
    97       CALL wrk_alloc( jpi, jpj, zswitch ) 
    98       CALL wrk_alloc( jpi, jpj, zhemis ) 
    99       CALL wrk_alloc( jpl,   2, zh_i_ini,  za_i_ini,  zv_i_ini ) 
    100       CALL wrk_alloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 
     102      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv  
     103      REAL(wp), POINTER, DIMENSION(:,:)   :: zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file 
     104      REAL(wp), POINTER, DIMENSION(:,:)   :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 
     105      REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini, zv_i_ini               !data by cattegories to fill 
     106      !-------------------------------------------------------------------- 
     107 
     108      CALL wrk_alloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini,  zv_i_ini ) 
     109      CALL wrk_alloc( jpi, jpj,      zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 
    101110 
    102111      IF(lwp) WRITE(numout,*) 
     
    117126 
    118127      ! basal temperature (considered at freezing point) 
    119       CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 
    120       t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)  
     128      t_bo(:,:) = ( eos_fzp( sss_m(:,:) ) + rt0 ) * tmask(:,:,1)  
    121129 
    122130      IF( ln_iceini ) THEN 
    123131 
    124       !-------------------------------------------------------------------- 
    125       ! 2) Basal temperature, ice mask and hemispheric index 
    126       !-------------------------------------------------------------------- 
    127  
    128       DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest 
    129          DO ji = 1, jpi 
    130             IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN  
    131                zswitch(ji,jj) = 0._wp * tmask(ji,jj,1)    ! no ice 
    132             ELSE                                                                                    
    133                zswitch(ji,jj) = 1._wp * tmask(ji,jj,1)    !    ice 
    134             ENDIF 
     132         !-------------------------------------------------------------------- 
     133         ! 2) Basal temperature, ice mask and hemispheric index 
     134         !-------------------------------------------------------------------- 
     135 
     136         DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest 
     137            DO ji = 1, jpi 
     138               IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN  
     139                  zswitch(ji,jj) = 0._wp * tmask(ji,jj,1)    ! no ice 
     140               ELSE                                                                                    
     141                  zswitch(ji,jj) = 1._wp * tmask(ji,jj,1)    !    ice 
     142               ENDIF 
     143            END DO 
    135144         END DO 
    136       END DO 
    137  
    138  
    139       ! Hemispheric index 
    140       DO jj = 1, jpj 
    141          DO ji = 1, jpi 
    142             IF( fcor(ji,jj) >= 0._wp ) THEN     
    143                zhemis(ji,jj) = 1 ! Northern hemisphere 
    144             ELSE 
    145                zhemis(ji,jj) = 2 ! Southern hemisphere 
    146             ENDIF 
    147          END DO 
    148       END DO 
    149  
    150       !-------------------------------------------------------------------- 
    151       ! 3) Initialization of sea ice state variables 
    152       !-------------------------------------------------------------------- 
    153  
    154       !----------------------------- 
    155       ! 3.1) Hemisphere-dependent arrays 
    156       !----------------------------- 
    157       ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 
    158       zht_i_ini(1) = rn_hti_ini_n ; zht_i_ini(2) = rn_hti_ini_s  ! ice thickness 
    159       zht_s_ini(1) = rn_hts_ini_n ; zht_s_ini(2) = rn_hts_ini_s  ! snow depth 
    160       zat_i_ini(1) = rn_ati_ini_n ; zat_i_ini(2) = rn_ati_ini_s  ! ice concentration 
    161       zsm_i_ini(1) = rn_smi_ini_n ; zsm_i_ini(2) = rn_smi_ini_s  ! bulk ice salinity 
    162       ztm_i_ini(1) = rn_tmi_ini_n ; ztm_i_ini(2) = rn_tmi_ini_s  ! temperature (ice and snow) 
    163  
    164       zvt_i_ini(:) = zht_i_ini(:) * zat_i_ini(:)   ! ice volume 
    165  
    166       !--------------------------------------------------------------------- 
    167       ! 3.2) Distribute ice concentration and thickness into the categories 
    168       !--------------------------------------------------------------------- 
    169       ! a gaussian distribution for ice concentration is used 
    170       ! then we check whether the distribution fullfills 
    171       ! volume and area conservation, positivity and ice categories bounds 
    172       DO i_hemis = 1, 2 
    173  
    174       ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 
    175  
    176       ! note for the great nemo engineers:  
    177       ! only very few of the WRITE statements are necessary for the reference version 
    178       ! they were one day useful, but now i personally doubt of their 
    179       ! potential for bringing anything useful 
    180  
    181       DO i_fill = jpl, 1, -1 
    182          IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN 
    183             !---------------------------- 
    184             ! fill the i_fill categories 
    185             !---------------------------- 
    186             ! *** 1 category to fill 
    187             IF ( i_fill .EQ. 1 ) THEN 
    188                zh_i_ini(1,i_hemis)       = zht_i_ini(i_hemis) 
    189                za_i_ini(1,i_hemis)       = zat_i_ini(i_hemis) 
    190                zh_i_ini(2:jpl,i_hemis)   = 0._wp 
    191                za_i_ini(2:jpl,i_hemis)   = 0._wp 
    192             ELSE 
    193  
    194                ! *** >1 categores to fill 
    195                !--- Ice thicknesses in the i_fill - 1 first categories 
    196                DO jl = 1, i_fill - 1 
    197                   zh_i_ini(jl,i_hemis) = hi_mean(jl) 
    198                END DO 
    199                 
    200                !--- jl0: most likely index where cc will be maximum 
    201                DO jl = 1, jpl 
    202                   IF ( ( zht_i_ini(i_hemis) >  hi_max(jl-1) ) .AND. & 
    203                      & ( zht_i_ini(i_hemis) <= hi_max(jl)   ) ) THEN 
    204                      jl0 = jl 
     145 
     146         !-------------------------------------------------------------------- 
     147         ! 3) Initialization of sea ice state variables 
     148         !-------------------------------------------------------------------- 
     149         IF( ln_iceini_file )THEN 
     150 
     151            zht_i_ini(:,:)  = si(jp_hti)%fnow(:,:,1) 
     152            zht_s_ini(:,:)  = si(jp_hts)%fnow(:,:,1) 
     153            zat_i_ini(:,:)  = si(jp_ati)%fnow(:,:,1) 
     154            zts_u_ini(:,:)  = si(jp_tsu)%fnow(:,:,1) 
     155            ztm_i_ini(:,:)  = si(jp_tmi)%fnow(:,:,1) 
     156            zsm_i_ini(:,:)  = si(jp_smi)%fnow(:,:,1) 
     157 
     158         ELSE ! ln_iceini_file = F 
     159 
     160            !----------------------------- 
     161            ! 3.1) Hemisphere-dependent arrays 
     162            !----------------------------- 
     163            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 
     164            DO jj = 1, jpj 
     165               DO ji = 1, jpi 
     166                  IF( fcor(ji,jj) >= 0._wp ) THEN 
     167                     zht_i_ini(ji,jj) = rn_hti_ini_n 
     168                     zht_s_ini(ji,jj) = rn_hts_ini_n 
     169                     zat_i_ini(ji,jj) = rn_ati_ini_n 
     170                     zsm_i_ini(ji,jj) = rn_smi_ini_n 
     171                     ztm_i_ini(ji,jj) = rn_tmi_ini_n 
     172                  ELSE 
     173                     zht_i_ini(ji,jj) = rn_hti_ini_s 
     174                     zht_s_ini(ji,jj) = rn_hts_ini_s 
     175                     zat_i_ini(ji,jj) = rn_ati_ini_s 
     176                     zsm_i_ini(ji,jj) = rn_smi_ini_s 
     177                     ztm_i_ini(ji,jj) = rn_tmi_ini_s 
    205178                  ENDIF 
    206179               END DO 
    207                jl0 = MIN(jl0, i_fill) 
    208                 
    209                !--- Concentrations 
    210                za_i_ini(jl0,i_hemis)      = zat_i_ini(i_hemis) / SQRT(REAL(jpl)) 
    211                DO jl = 1, i_fill - 1 
    212                   IF ( jl .NE. jl0 ) THEN 
    213                      zsigma               = 0.5 * zht_i_ini(i_hemis) 
    214                      zarg                 = ( zh_i_ini(jl,i_hemis) - zht_i_ini(i_hemis) ) / zsigma 
    215                      za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2) 
    216                   ENDIF 
    217                END DO 
    218                 
    219                zA = 0. ! sum of the areas in the jpl categories  
    220                DO jl = 1, i_fill - 1 
    221                  zA = zA + za_i_ini(jl,i_hemis) 
    222                END DO 
    223                za_i_ini(i_fill,i_hemis)   = zat_i_ini(i_hemis) - zA ! ice conc in the last category 
    224                IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 
    225           
    226                !--- Ice thickness in the last category 
    227                zV = 0. ! sum of the volumes of the N-1 categories 
    228                DO jl = 1, i_fill - 1 
    229                   zV = zV + za_i_ini(jl,i_hemis)*zh_i_ini(jl,i_hemis) 
    230                END DO 
    231                zh_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis)  
    232                IF ( i_fill .LT. jpl ) zh_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 
    233  
    234                !--- volumes 
    235                zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zh_i_ini(:,i_hemis) 
    236                IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 
    237  
    238             ENDIF ! i_fill 
    239  
    240             !--------------------- 
    241             ! Compatibility tests 
    242             !--------------------- 
    243             ! Test 1: area conservation 
    244             zA_cons = SUM(za_i_ini(:,i_hemis)) ; zconv = ABS(zat_i_ini(i_hemis) - zA_cons ) 
    245             IF ( zconv .LT. 1.0e-6 ) THEN 
    246                ztest_1 = 1 
    247             ELSE  
    248               ! this write is useful 
    249               IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis)  
    250                ztest_1 = 0 
    251             ENDIF 
    252  
    253             ! Test 2: volume conservation 
    254             zV_cons = SUM(zv_i_ini(:,i_hemis)) 
    255             zconv = ABS(zvt_i_ini(i_hemis) - zV_cons) 
    256  
    257             IF ( zconv .LT. 1.0e-6 ) THEN 
    258                ztest_2 = 1 
    259             ELSE 
    260               ! this write is useful 
    261               IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, & 
    262                             ' zvt_i_ini = ', zvt_i_ini(i_hemis) 
    263                ztest_2 = 0 
    264             ENDIF 
    265  
    266             ! Test 3: thickness of the last category is in-bounds ? 
    267             IF ( zh_i_ini(i_fill, i_hemis) > hi_max(i_fill-1) ) THEN 
    268                ztest_3 = 1 
    269             ELSE 
    270                ! this write is useful 
    271                IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(i_fill,i_hemis) = ', & 
    272                zh_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 
    273                ztest_3 = 0 
    274             ENDIF 
    275  
    276             ! Test 4: positivity of ice concentrations 
    277             ztest_4 = 1 
    278             DO jl = 1, jpl 
    279                IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN  
    280                   ! this write is useful 
    281                   IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis) 
    282                   ztest_4 = 0 
    283                ENDIF 
    284             END DO 
    285  
    286          ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 
    287   
    288          ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 
    289  
    290       END DO ! i_fill 
    291  
    292       IF(lwp) THEN  
    293          WRITE(numout,*) ' ztests : ', ztests 
    294          IF ( ztests .NE. 4 ) THEN 
    295             WRITE(numout,*) 
    296             WRITE(numout,*) ' !!!! ALERT                  !!! ' 
    297             WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure ' 
    298             WRITE(numout,*) 
    299             WRITE(numout,*) ' *** ztests is not equal to 4 ' 
    300             WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 
    301             WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(i_hemis) 
    302             WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(i_hemis) 
    303          ENDIF ! ztests .NE. 4 
    304       ENDIF 
    305        
    306       END DO ! i_hemis 
    307  
    308       !--------------------------------------------------------------------- 
    309       ! 3.3) Space-dependent arrays for ice state variables 
    310       !--------------------------------------------------------------------- 
    311  
    312       ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature 
    313       DO jl = 1, jpl ! loop over categories 
     180            END DO 
     181 
     182         ENDIF ! ln_iceini_file 
     183 
     184         zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:)   ! ice volume 
     185 
     186         !--------------------------------------------------------------------- 
     187         ! 3.2) Distribute ice concentration and thickness into the categories 
     188         !--------------------------------------------------------------------- 
     189         ! a gaussian distribution for ice concentration is used 
     190         ! then we check whether the distribution fullfills 
     191         ! volume and area conservation, positivity and ice categories bounds 
     192         zh_i_ini(:,:,:) = 0._wp  
     193         za_i_ini(:,:,:) = 0._wp 
     194         zv_i_ini(:,:,:) = 0._wp 
     195 
    314196         DO jj = 1, jpj 
    315197            DO ji = 1, jpi 
    316                a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration 
    317                ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj))   ! ice thickness 
    318                ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zht_s_ini( zhemis(ji,jj) ) / zht_i_ini( zhemis(ji,jj) ) )  ! snow depth 
    319                sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj))     ! salinity 
    320                o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp                        ! age (1 day) 
    321                t_su(ji,jj,jl)  = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp 
    322  
    323                ! This case below should not be used if (ht_s/ht_i) is ok in namelist 
    324                ! In case snow load is in excess that would lead to transformation from snow to ice 
    325                ! Then, transfer the snow excess into the ice (different from limthd_dh) 
    326                zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 )  
    327                ! recompute ht_i, ht_s avoiding out of bounds values 
    328                ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 
    329                ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn ) 
    330  
    331                ! ice volume, salt content, age content 
    332                v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume 
    333                v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume 
    334                smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 
    335                oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content 
    336             END DO 
    337          END DO 
    338       END DO 
    339  
    340       ! Snow temperature and heat content 
    341       DO jk = 1, nlay_s 
     198 
     199               IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN 
     200 
     201                  ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 
     202                  ztests  = 0  
     203 
     204                  DO i_fill = jpl, 1, -1 
     205 
     206                     IF( ztests .NE. 4 ) THEN 
     207                        !---------------------------- 
     208                        ! fill the i_fill categories 
     209                        !---------------------------- 
     210                        ! *** 1 category to fill 
     211                        IF ( i_fill .EQ. 1 ) THEN 
     212                           zh_i_ini(ji,jj,    1)   = zht_i_ini(ji,jj) 
     213                           za_i_ini(ji,jj,    1)   = zat_i_ini(ji,jj) 
     214                           zh_i_ini(ji,jj,2:jpl)   = 0._wp 
     215                           za_i_ini(ji,jj,2:jpl)   = 0._wp 
     216                        ELSE 
     217 
     218                           ! *** >1 categores to fill 
     219                           !--- Ice thicknesses in the i_fill - 1 first categories 
     220                           DO jl = 1, i_fill - 1 
     221                              zh_i_ini(ji,jj,jl) = hi_mean(jl) 
     222                           END DO 
     223                
     224                           !--- jl0: most likely index where cc will be maximum 
     225                           DO jl = 1, jpl 
     226                              IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. & 
     227                                 & ( zht_i_ini(ji,jj) <= hi_max(jl)   ) ) THEN 
     228                                 jl0 = jl 
     229                              ENDIF 
     230                           END DO 
     231                           jl0 = MIN(jl0, i_fill) 
     232                
     233                           !--- Concentrations 
     234                           za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 
     235                           DO jl = 1, i_fill - 1 
     236                              IF( jl .NE. jl0 )THEN 
     237                                 zsigma             = 0.5 * zht_i_ini(ji,jj) 
     238                                 zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / zsigma 
     239                                 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 
     240                              ENDIF 
     241                           END DO 
     242                
     243                           zA = 0. ! sum of the areas in the jpl categories  
     244                           DO jl = 1, i_fill - 1 
     245                              zA = zA + za_i_ini(ji,jj,jl) 
     246                           END DO 
     247                           za_i_ini(ji,jj,i_fill)   = zat_i_ini(ji,jj) - zA ! ice conc in the last category 
     248                           IF ( i_fill .LT. jpl ) za_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 
     249          
     250                           !--- Ice thickness in the last category 
     251                           zV = 0. ! sum of the volumes of the N-1 categories 
     252                           DO jl = 1, i_fill - 1 
     253                              zV = zV + za_i_ini(ji,jj,jl)*zh_i_ini(ji,jj,jl) 
     254                           END DO 
     255                           zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / za_i_ini(ji,jj,i_fill)  
     256                           IF ( i_fill .LT. jpl ) zh_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 
     257 
     258                           !--- volumes 
     259                           zv_i_ini(ji,jj,:) = za_i_ini(ji,jj,:) * zh_i_ini(ji,jj,:) 
     260                           IF ( i_fill .LT. jpl ) zv_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 
     261 
     262                        ENDIF ! i_fill 
     263 
     264                        !--------------------- 
     265                        ! Compatibility tests 
     266                        !--------------------- 
     267                        ! Test 1: area conservation 
     268                        zA_cons = SUM(za_i_ini(ji,jj,:)) ; zconv = ABS(zat_i_ini(ji,jj) - zA_cons ) 
     269                        IF ( zconv .LT. 1.0e-6 ) THEN 
     270                           ztest_1 = 1 
     271                        ELSE  
     272                          !this write is useful 
     273                          IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(ji,jj)  
     274                          ztest_1 = 0 
     275                        ENDIF 
     276 
     277                        ! Test 2: volume conservation 
     278                        zV_cons = SUM(zv_i_ini(ji,jj,:)) 
     279                        zconv = ABS(zvt_i_ini(ji,jj) - zV_cons) 
     280 
     281                        IF( zconv .LT. 1.0e-6 ) THEN 
     282                           ztest_2 = 1 
     283                        ELSE 
     284                           !this write is useful 
     285                           IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, & 
     286                                                    ' zvt_i_ini = ', zvt_i_ini(ji,jj) 
     287                           ztest_2 = 0 
     288                        ENDIF 
     289 
     290                        ! Test 3: thickness of the last category is in-bounds ? 
     291                        IF ( zh_i_ini(ji,jj,i_fill) > hi_max(i_fill-1) ) THEN 
     292                           ztest_3 = 1 
     293                        ELSE 
     294                           ! this write is useful 
     295                           IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(ji,jj,i_fill) = ', & 
     296                           zh_i_ini(ji,jj,i_fill), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 
     297                           ztest_3 = 0 
     298                        ENDIF 
     299 
     300                        ! Test 4: positivity of ice concentrations 
     301                        ztest_4 = 1 
     302                        DO jl = 1, jpl 
     303                           IF ( za_i_ini(ji,jj,jl) .LT. 0._wp ) THEN  
     304                              ! this write is useful 
     305                              IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(ji,jj,jl) 
     306                              ztest_4 = 0 
     307                           ENDIF 
     308                        END DO 
     309 
     310                     ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 
     311  
     312                     ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 
     313 
     314                  END DO ! i_fill 
     315 
     316                  IF(lwp) THEN  
     317                     WRITE(numout,*) ' ztests : ', ztests 
     318                     IF( ztests .NE. 4 )THEN 
     319                        WRITE(numout,*) 
     320                        WRITE(numout,*) ' !!!! ALERT                  !!! ' 
     321                        WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure ' 
     322                        WRITE(numout,*) 
     323                        WRITE(numout,*) ' *** ztests is not equal to 4 ' 
     324                        WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 
     325                        WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 
     326                        WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 
     327                     ENDIF ! ztests .NE. 4 
     328                  ENDIF 
     329       
     330               ENDIF !  zat_i_ini(ji,jj) > 0._wp .AND. zhm_i_ini(ji,jj) > 0._wp 
     331 
     332            ENDDO    
     333         ENDDO    
     334 
     335         !--------------------------------------------------------------------- 
     336         ! 3.3) Space-dependent arrays for ice state variables 
     337         !--------------------------------------------------------------------- 
     338 
     339         ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature 
    342340         DO jl = 1, jpl ! loop over categories 
    343341            DO jj = 1, jpj 
    344342               DO ji = 1, jpi 
    345                    t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 
    346                    ! Snow energy of melting 
    347                    e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 
    348  
    349                    ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 
    350                    e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 
     343                  a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini(ji,jj,jl)                       ! concentration 
     344                  ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(ji,jj,jl)                       ! ice thickness 
     345                  sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(ji,jj)                         ! salinity 
     346                  o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp                                    ! age (1 day) 
     347                  t_su(ji,jj,jl)  = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp 
     348 
     349                  IF( zht_i_ini(ji,jj) > 0._wp )THEN 
     350                    ht_s(ji,jj,jl)= ht_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) )  ! snow depth 
     351                  ELSE 
     352                    ht_s(ji,jj,jl)= 0._wp 
     353                  ENDIF 
     354 
     355                  ! This case below should not be used if (ht_s/ht_i) is ok in namelist 
     356                  ! In case snow load is in excess that would lead to transformation from snow to ice 
     357                  ! Then, transfer the snow excess into the ice (different from limthd_dh) 
     358                  zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 )  
     359                  ! recompute ht_i, ht_s avoiding out of bounds values 
     360                  ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 
     361                  ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn ) 
     362 
     363                  ! ice volume, salt content, age content 
     364                  v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume 
     365                  v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume 
     366                  smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 
     367                  oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content 
    351368               END DO 
    352369            END DO 
    353370         END DO 
    354       END DO 
    355  
    356       ! Ice salinity, temperature and heat content 
    357       DO jk = 1, nlay_i 
    358          DO jl = 1, jpl ! loop over categories 
    359             DO jj = 1, jpj 
    360                DO ji = 1, jpi 
    361                    t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0  
    362                    s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin 
    363                    ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 
    364  
    365                    ! heat content per unit volume 
    366                    e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
    367                       +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 
    368                       -   rcp     * ( ztmelts - rt0 ) ) 
    369  
    370                    ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
    371                    e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 
     371 
     372         ! Snow temperature and heat content 
     373         DO jk = 1, nlay_s 
     374            DO jl = 1, jpl ! loop over categories 
     375               DO jj = 1, jpj 
     376                  DO ji = 1, jpi 
     377                     t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 
     378                     ! Snow energy of melting 
     379                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 
     380 
     381                     ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 
     382                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 
     383                  END DO 
    372384               END DO 
    373385            END DO 
    374386         END DO 
    375       END DO 
    376  
    377       tn_ice (:,:,:) = t_su (:,:,:) 
    378  
    379       ELSE  
    380          ! if ln_iceini=false 
     387 
     388         ! Ice salinity, temperature and heat content 
     389         DO jk = 1, nlay_i 
     390            DO jl = 1, jpl ! loop over categories 
     391               DO jj = 1, jpj 
     392                  DO ji = 1, jpi 
     393                     t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0  
     394                     s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin 
     395                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 
     396 
     397                     ! heat content per unit volume 
     398                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
     399                        +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 
     400                        -   rcp     * ( ztmelts - rt0 ) ) 
     401 
     402                     ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
     403                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 
     404                  END DO 
     405               END DO 
     406            END DO 
     407         END DO 
     408 
     409         tn_ice (:,:,:) = t_su (:,:,:) 
     410 
     411      ELSE ! if ln_iceini=false 
    381412         a_i  (:,:,:) = 0._wp 
    382413         v_i  (:,:,:) = 0._wp 
     
    400431            END DO 
    401432         END DO 
    402        
     433 
    403434      ENDIF ! ln_iceini 
    404435       
     
    452483 
    453484 
    454       CALL wrk_dealloc( jpi, jpj, zswitch ) 
    455       CALL wrk_dealloc( jpi, jpj, zhemis ) 
    456       CALL wrk_dealloc( jpl,   2, zh_i_ini,  za_i_ini,  zv_i_ini ) 
    457       CALL wrk_dealloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 
     485      CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini,  zv_i_ini ) 
     486      CALL wrk_dealloc( jpi, jpj,      zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 
    458487 
    459488   END SUBROUTINE lim_istate 
     
    475504      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization 
    476505      !!----------------------------------------------------------------------------- 
    477       NAMELIST/namiceini/ ln_iceini, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, rn_hti_ini_n, rn_hti_ini_s,  & 
    478          &                                      rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s 
    479       INTEGER :: ios                 ! Local integer output status for namelist read 
     506      ! 
     507      INTEGER :: ios,ierr,inum_ice                 ! Local integer output status for namelist read 
     508      INTEGER :: ji,jj 
     509      INTEGER :: ifpr, ierror 
     510      ! 
     511      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ice files 
     512      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi 
     513      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read 
     514      ! 
     515      NAMELIST/namiceini/ ln_iceini, ln_iceini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s,  & 
     516         &                rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, & 
     517         &                rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s,                             & 
     518         &                sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir 
    480519      !!----------------------------------------------------------------------------- 
    481520      ! 
     
    489528      IF(lwm) WRITE ( numoni, namiceini ) 
    490529 
     530      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts 
     531      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu 
     532      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi 
     533 
    491534      ! Define the initial parameters 
    492535      ! ------------------------- 
     
    497540         WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    498541         WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_iceini     = ', ln_iceini 
     542         WRITE(numout,*) '   ice initialization from a netcdf file      ln_iceini_file  = ', ln_iceini_file 
    499543         WRITE(numout,*) '   threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst 
    500544         WRITE(numout,*) '   initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n 
     
    510554      ENDIF 
    511555 
     556      IF( ln_iceini_file ) THEN                      ! Ice initialization using input file 
     557         ! 
     558         ! set si structure 
     559         ALLOCATE( si(jpfldi), STAT=ierror ) 
     560         IF( ierror > 0 ) THEN 
     561            CALL ctl_stop( 'Ice_ini in limistate: unable to allocate si structure' )   ;   RETURN 
     562         ENDIF 
     563 
     564         DO ifpr = 1, jpfldi 
     565            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) ) 
     566            ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) ) 
     567         END DO 
     568 
     569         ! fill si with slf_i and control print 
     570         CALL fld_fill( si, slf_i, cn_dir, 'lim_istate', 'lim istate ini', 'numnam_ice' ) 
     571 
     572         CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step 
     573 
     574      ENDIF 
     575 
    512576   END SUBROUTINE lim_istate_init 
    513577 
Note: See TracChangeset for help on using the changeset viewer.