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

Changeset 6853


Ignore:
Timestamp:
2016-08-08T11:02:18+02:00 (8 years ago)
Author:
clem
Message:

correct bugs on LIM3 rheology and outputs. Make sure the ice thickness distribution (used in initialization and bdy boundary conditions) is a gaussian.

Location:
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r6746 r6853  
    3131   USE iom 
    3232 
     33!!!clem 
     34!!   USE diawri 
     35!!! 
     36    
    3337   IMPLICIT NONE 
    3438   PRIVATE 
     
    8488      REAL(wp)   :: ztmelts, zdh 
    8589      INTEGER    :: i_hemis, i_fill, jl0   
    86       REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv  
     90      REAL(wp)   :: zarg, zV, zconv, zdv  
    8791      REAL(wp), POINTER, DIMENSION(:,:)   :: zswitch    ! ice indicator 
    8892      REAL(wp), POINTER, DIMENSION(:,:)   :: zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file 
    8993      REAL(wp), POINTER, DIMENSION(:,:)   :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 
    90       REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini, zv_i_ini               !data by cattegories to fill 
    91       !-------------------------------------------------------------------- 
    92  
    93       CALL wrk_alloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini,  zv_i_ini ) 
     94      REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini                         !data by cattegories to fill 
     95      INTEGER , POINTER, DIMENSION(:)     :: itest 
     96      !-------------------------------------------------------------------- 
     97 
     98      CALL wrk_alloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini ) 
    9499      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 ) 
    95100      CALL wrk_alloc( jpi, jpj,      zswitch ) 
     101      Call wrk_alloc( 4,             itest ) 
    96102 
    97103      IF(lwp) WRITE(numout,*) 
     
    124130         DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest 
    125131            DO ji = 1, jpi 
    126                IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN  
    127                   zswitch(ji,jj) = 0._wp * tmask(ji,jj,1)    ! no ice 
     132               IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst * tmask(ji,jj,1) ) THEN 
     133                  zswitch(ji,jj) = 0._wp                     ! no ice 
    128134               ELSE                                                                                    
    129135                  zswitch(ji,jj) = 1._wp * tmask(ji,jj,1)    !    ice 
     
    153159               DO ji = 1, jpi 
    154160                  IF( ff(ji,jj) >= 0._wp ) THEN 
    155                      zht_i_ini(ji,jj) = rn_hti_ini_n 
    156                      zht_s_ini(ji,jj) = rn_hts_ini_n 
    157                      zat_i_ini(ji,jj) = rn_ati_ini_n 
    158                      zts_u_ini(ji,jj) = rn_tmi_ini_n 
    159                      zsm_i_ini(ji,jj) = rn_smi_ini_n 
    160                      ztm_i_ini(ji,jj) = rn_tmi_ini_n 
     161                     zht_i_ini(ji,jj) = rn_hti_ini_n * zswitch(ji,jj) 
     162                     zht_s_ini(ji,jj) = rn_hts_ini_n * zswitch(ji,jj) 
     163                     zat_i_ini(ji,jj) = rn_ati_ini_n * zswitch(ji,jj) 
     164                     zts_u_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj) 
     165                     zsm_i_ini(ji,jj) = rn_smi_ini_n * zswitch(ji,jj) 
     166                     ztm_i_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj) 
    161167                  ELSE 
    162                      zht_i_ini(ji,jj) = rn_hti_ini_s 
    163                      zht_s_ini(ji,jj) = rn_hts_ini_s 
    164                      zat_i_ini(ji,jj) = rn_ati_ini_s 
    165                      zts_u_ini(ji,jj) = rn_tmi_ini_s 
    166                      zsm_i_ini(ji,jj) = rn_smi_ini_s 
    167                      ztm_i_ini(ji,jj) = rn_tmi_ini_s 
     168                     zht_i_ini(ji,jj) = rn_hti_ini_s * zswitch(ji,jj) 
     169                     zht_s_ini(ji,jj) = rn_hts_ini_s * zswitch(ji,jj) 
     170                     zat_i_ini(ji,jj) = rn_ati_ini_s * zswitch(ji,jj) 
     171                     zts_u_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj) 
     172                     zsm_i_ini(ji,jj) = rn_smi_ini_s * zswitch(ji,jj) 
     173                     ztm_i_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj) 
    168174                  ENDIF 
    169175               END DO 
     
    181187         zh_i_ini(:,:,:) = 0._wp  
    182188         za_i_ini(:,:,:) = 0._wp 
    183          zv_i_ini(:,:,:) = 0._wp 
    184189 
    185190         DO jj = 1, jpj 
     
    188193               IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN 
    189194 
    190                   ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 
    191 !                  ztests  = 0  
    192  
    193                   DO i_fill = jpl, 1, -1 
    194  
    195 !                     IF( ztests .NE. 4 ) THEN 
    196                      IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN 
    197                         !---------------------------- 
    198                         ! fill the i_fill categories 
    199                         !---------------------------- 
    200                         ! *** 1 category to fill 
    201                         IF ( i_fill .EQ. 1 ) THEN 
    202                            zh_i_ini(ji,jj,    1)   = zht_i_ini(ji,jj) 
    203                            za_i_ini(ji,jj,    1)   = zat_i_ini(ji,jj) 
    204                            zh_i_ini(ji,jj,2:jpl)   = 0._wp 
    205                            za_i_ini(ji,jj,2:jpl)   = 0._wp 
    206                         ELSE 
    207  
    208                            ! *** >1 categores to fill 
    209                            !--- Ice thicknesses in the i_fill - 1 first categories 
    210                            DO jl = 1, i_fill - 1 
    211                               zh_i_ini(ji,jj,jl) = hi_mean(jl) 
    212                            END DO 
     195                  !--- jl0: most likely index where cc will be maximum 
     196                  jl0 = jpl 
     197                  DO jl = 1, jpl 
     198                     IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN 
     199                        jl0 = jl 
     200                        CYCLE 
     201                     ENDIF 
     202                  END DO 
     203 
     204                  ! initialisation of tests 
     205                  itest(:)  = 0 
     206                   
     207                  i_fill = jpl + 1                                             !==================================== 
     208                  DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories 
     209                     ! iteration                                               !==================================== 
     210                     i_fill = i_fill - 1 
     211 
     212                     ! initialisation of ice variables for each try 
     213                     zh_i_ini(ji,jj,:) = 0._wp  
     214                     za_i_ini(ji,jj,:) = 0._wp 
     215                     itest(:) = 0 
     216 
     217                     ! *** case very thin ice: fill only category 1 
     218                     IF ( i_fill == 1 ) THEN 
     219                        zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj) 
     220                        za_i_ini(ji,jj,1) = zat_i_ini(ji,jj) 
     221 
     222                     ! *** case ice is thicker: fill categories >1 
     223                     ELSE 
     224 
     225                        ! Fill ice thicknesses in the (i_fill-1) cat by hmean  
     226                        DO jl = 1, i_fill-1 
     227                           zh_i_ini(ji,jj,jl) = hi_mean(jl) 
     228                        END DO 
    213229                
    214                            !--- jl0: most likely index where cc will be maximum 
    215                            DO jl = 1, jpl 
    216                               IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. & 
    217                                  & ( zht_i_ini(ji,jj) <= hi_max(jl)   ) ) THEN 
    218                                  jl0 = jl 
    219                               ENDIF 
    220                            END DO 
    221                            jl0 = MIN(jl0, i_fill) 
    222                 
    223                            !--- Concentrations 
    224                            za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 
    225                            DO jl = 1, i_fill - 1 
    226                               IF( jl .NE. jl0 )THEN 
    227                                  zsigma             = 0.5 * zht_i_ini(ji,jj) 
    228                                  zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / zsigma 
    229                                  za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 
    230                               ENDIF 
    231                            END DO 
    232                 
    233                            zA = 0. ! sum of the areas in the jpl categories  
    234                            DO jl = 1, i_fill - 1 
    235                               zA = zA + za_i_ini(ji,jj,jl) 
    236                            END DO 
    237                            za_i_ini(ji,jj,i_fill)   = zat_i_ini(ji,jj) - zA ! ice conc in the last category 
    238                            IF ( i_fill .LT. jpl ) za_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 
    239           
    240                            !--- Ice thickness in the last category 
    241                            zV = 0. ! sum of the volumes of the N-1 categories 
    242                            DO jl = 1, i_fill - 1 
    243                               zV = zV + za_i_ini(ji,jj,jl)*zh_i_ini(ji,jj,jl) 
    244                            END DO 
    245                            zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / za_i_ini(ji,jj,i_fill)  
    246                            IF ( i_fill .LT. jpl ) zh_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 
    247  
    248                            !--- volumes 
    249                            zv_i_ini(ji,jj,:) = za_i_ini(ji,jj,:) * zh_i_ini(ji,jj,:) 
    250                            IF ( i_fill .LT. jpl ) zv_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 
    251  
    252                         ENDIF ! i_fill 
    253  
    254                         !--------------------- 
    255                         ! Compatibility tests 
    256                         !--------------------- 
    257                         ! Test 1: area conservation 
    258                         zA_cons = SUM(za_i_ini(ji,jj,:)) ; zconv = ABS(zat_i_ini(ji,jj) - zA_cons ) 
    259                         IF ( zconv .LT. 1.0e-6 ) THEN 
    260                            ztest_1 = 1 
    261                         ELSE  
    262                           ztest_1 = 0 
    263                         ENDIF 
    264  
    265                         ! Test 2: volume conservation 
    266                         zV_cons = SUM(zv_i_ini(ji,jj,:)) 
    267                         zconv = ABS(zvt_i_ini(ji,jj) - zV_cons) 
    268  
    269                         IF( zconv .LT. 1.0e-6 ) THEN 
    270                            ztest_2 = 1 
    271                         ELSE 
    272                            ztest_2 = 0 
    273                         ENDIF 
    274  
    275                         ! Test 3: thickness of the last category is in-bounds ? 
    276                         IF ( zh_i_ini(ji,jj,i_fill) > hi_max(i_fill-1) ) THEN 
    277                            ztest_3 = 1 
    278                         ELSE 
    279                            ztest_3 = 0 
    280                         ENDIF 
    281  
    282                         ! Test 4: positivity of ice concentrations 
    283                         ztest_4 = 1 
    284                         DO jl = 1, jpl 
    285                            IF ( za_i_ini(ji,jj,jl) .LT. 0._wp ) THEN  
    286                               ztest_4 = 0 
     230                        !--- Concentrations 
     231                        za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 
     232                        DO jl = 1, i_fill - 1 
     233                           IF( jl /= jl0 )THEN 
     234                              zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) ) 
     235                              za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 
    287236                           ENDIF 
    288237                        END DO 
    289  
    290                      ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 
    291   
    292                      ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 
    293  
    294                   END DO ! i_fill 
    295  
    296                   IF(lwp) THEN  
    297                      WRITE(numout,*) ' ztests : ', ztests 
    298                      IF( ztests .NE. 4 )THEN 
    299                         WRITE(numout,*) 
    300                         WRITE(numout,*) ' !!!! ALERT                  !!! ' 
    301                         WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure ' 
    302                         WRITE(numout,*) 
    303                         WRITE(numout,*) ' *** ztests is not equal to 4 ' 
    304                         WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 
    305                         WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 
    306                         WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 
    307                      ENDIF ! ztests .NE. 4 
     238                
     239                        ! Concentration in the last (i_fill) category 
     240                        za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) ) 
     241 
     242                        ! Ice thickness in the last (i_fill) category 
     243                        zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) ) 
     244                        zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 )  
     245 
     246                        ! clem: correction if concentration of upper cat is greater than lower cat 
     247                        !       (it should be a gaussian around jl0 but sometimes it is not) 
     248                        IF ( jl0 /= jpl ) THEN 
     249                           DO jl = jpl, jl0+1, -1 
     250                              IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN 
     251                                 zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl) 
     252                                 zh_i_ini(ji,jj,jl    ) = 0._wp 
     253                                 za_i_ini(ji,jj,jl    ) = 0._wp 
     254                                 za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  & 
     255                                    &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 ) 
     256                              END IF 
     257                           ENDDO 
     258                        ENDIF 
     259 
     260                     ENDIF ! case ice is thick or thin 
     261 
     262                     !--------------------- 
     263                     ! Compatibility tests 
     264                     !--------------------- 
     265                     ! Test 1: area conservation 
     266                     zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) ) 
     267                     IF ( zconv < epsi06 ) itest(1) = 1 
     268                      
     269                     ! Test 2: volume conservation 
     270                     zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   & 
     271                        &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) ) 
     272                     IF ( zconv < epsi06 ) itest(2) = 1 
     273                      
     274                     ! Test 3: thickness of the last category is in-bounds ? 
     275                     IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 
     276                      
     277                     ! Test 4: positivity of ice concentrations 
     278                     itest(4) = 1 
     279                     DO jl = 1, i_fill 
     280                        IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0 
     281                     END DO 
     282                     !                                      !============================ 
     283                  END DO                                    ! end iteration on categories 
     284                  !                                         !============================ 
     285 
     286                  IF( lwp .AND. SUM(itest) /= 4 ) THEN  
     287                     WRITE(numout,*) 
     288                     WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! ' 
     289                     WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure ' 
     290                     WRITE(numout,*) 
     291                     WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:) 
     292                     WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 
     293                     WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 
    308294                  ENDIF 
    309        
    310                ENDIF !  zat_i_ini(ji,jj) > 0._wp .AND. zhm_i_ini(ji,jj) > 0._wp 
    311  
    312             ENDDO    
    313          ENDDO    
     295                
     296               ENDIF !  zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp 
     297             
     298            ENDDO 
     299         ENDDO 
    314300 
    315301         !--------------------------------------------------------------------- 
     
    468454      sxyage (:,:,:)  = 0._wp 
    469455 
    470  
    471       CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini,  zv_i_ini ) 
     456!!!clem 
     457!!      ! Output the initial state and forcings 
     458!!      CALL dia_wri_state( 'output.init', nit000 ) 
     459!!! 
     460       
     461      CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini ) 
    472462      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 ) 
    473463      CALL wrk_dealloc( jpi, jpj,      zswitch ) 
     464      Call wrk_dealloc( 4,             itest ) 
    474465 
    475466   END SUBROUTINE lim_istate 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r6801 r6853  
    3838   USE bdyice_lim 
    3939#endif 
     40   !!!clem 
     41!!   USE iom 
    4042 
    4143   IMPLICIT NONE 
     
    140142                                                                             !   ocean surface (ssh_m) if ice is not embedded 
    141143                                                                             !   ice top surface if ice is embedded    
    142       REAL(wp), POINTER, DIMENSION(:,:) ::   zswitch1, zswitch2              ! dummy arrays 
     144      REAL(wp), POINTER, DIMENSION(:,:) ::   zswitchU, zswitchV              ! dummy arrays 
     145      REAL(wp), POINTER, DIMENSION(:,:) ::   zmaskU, zmaskV                  ! mask for ice presence 
    143146      REAL(wp), POINTER, DIMENSION(:,:) ::   zfmask, zwf                     ! mask at F points for the ice 
    144147 
    145148      REAL(wp), PARAMETER               ::   zepsi  = 1.0e-20_wp             ! tolerance parameter 
    146       REAL(wp), PARAMETER               ::   zvmin  = 1.0e-03_wp             ! ice volume below which ice velocity equals ocean velocity 
     149      REAL(wp), PARAMETER               ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity 
    147150      REAL(wp), PARAMETER               ::   zshlat = 2._wp                  ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip) 
    148151      !!------------------------------------------------------------------- 
     
    152155      CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV ) 
    153156      CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
    154       CALL wrk_alloc( jpi,jpj, zswitch1, zswitch2, zfmask, zwf ) 
     157      CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf ) 
    155158 
    156159#if defined key_agrif  
     
    294297            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 
    295298 
     299            ! masks 
     300            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice 
     301            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice 
     302 
    296303            ! switches 
    297             zswitch1(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) 
    298             zswitch2(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) 
     304            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin 
     305            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin 
    299306 
    300307         END DO 
     
    436443                   
    437444                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
    438                   v_ice(ji,jj) = (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity 
    439                      &                                   + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    440                      &                                   ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
    441                      &           + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    442                      &           ) * zswitch2(ji,jj) 
     445                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity 
     446                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     447                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
     448                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     449                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
     450                     &           ) * zmaskV(ji,jj) 
    443451                  ! Bouillon 2013 
    444452                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
    445453                  !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
    446                   !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 
     454                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 
    447455 
    448456               END DO 
     
    480488 
    481489                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
    482                   u_ice(ji,jj) = (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity 
    483                      &                                   + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    484                      &                                   ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
    485                      &           + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    486                      &           ) * zswitch1(ji,jj) 
     490                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity 
     491                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     492                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
     493                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     494                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin  
     495                     &           ) * zmaskU(ji,jj) 
    487496                  ! Bouillon 2013 
    488497                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
    489498                  !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
    490                   !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 
     499                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 
    491500               END DO 
    492501            END DO 
     
    525534 
    526535                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
    527                   u_ice(ji,jj) = (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity 
    528                      &                                   + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    529                      &                                   ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )            &  ! m/dt + tau_io(only ice part) + landfast 
    530                      &           + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    531                      &           ) * zswitch1(ji,jj) 
     536                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity 
     537                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     538                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
     539                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     540                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
     541                     &           ) * zmaskU(ji,jj) 
    532542                  ! Bouillon 2013 
    533543                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
    534544                  !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
    535                   !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 
     545                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 
    536546               END DO 
    537547            END DO 
     
    564574                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
    565575 
    566                   ! landfast switch => 0 = static friction ; 1 = sliding friction 
    567                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     576                  ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction 
     577                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    568578                   
    569579                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
    570                   v_ice(ji,jj) = (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity 
    571                      &                                   + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    572                      &                                   ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )            &  ! m/dt + tau_io(only ice part) + landfast 
    573                      &           + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    574                      &           ) * zswitch2(ji,jj) 
     580                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity 
     581                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     582                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
     583                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     584                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
     585                     &           ) * zmaskV(ji,jj) 
    575586                  ! Bouillon 2013 
    576587                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
    577588                  !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
    578                   !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 
     589                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 
    579590               END DO 
    580591            END DO 
     
    603614      ! 
    604615      !------------------------------------------------------------------------------! 
    605       ! 4) Limit ice velocities when ice is thin 
    606       !------------------------------------------------------------------------------! 
    607       ! If the ice volume is below zvmin then ice velocity should equal the 
    608       ! ocean velocity. This prevents high velocity when ice is thin 
    609        DO jj = 2, jpjm1 
    610          DO ji = fs_2, fs_jpim1 
    611             rswitch      = MAX( 0._wp, SIGN( 1._wp, vt_i(ji,jj) - zvmin ) ) 
    612             u_ice(ji,jj) = ( rswitch * u_ice(ji,jj) + (1._wp - rswitch) * u_oce(ji,jj) ) * zswitch1(ji,jj) 
    613             v_ice(ji,jj) = ( rswitch * v_ice(ji,jj) + (1._wp - rswitch) * v_oce(ji,jj) ) * zswitch2(ji,jj) 
    614          END DO 
    615       END DO 
    616  
    617       CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. ) 
    618  
    619 #if defined key_agrif 
    620       CALL agrif_interp_lim3( 'U' ) 
    621       CALL agrif_interp_lim3( 'V' ) 
    622 #endif 
    623 #if defined key_bdy 
    624       CALL bdy_ice_lim_dyn( 'U' ) 
    625       CALL bdy_ice_lim_dyn( 'V' ) 
    626 #endif          
    627  
    628       !------------------------------------------------------------------------------! 
    629       ! 5) Recompute delta, shear and div (inputs for mechanical redistribution)  
     616      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)  
    630617      !------------------------------------------------------------------------------! 
    631618      DO jj = 1, jpjm1 
     
    665652            ! delta at T points 
    666653            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 )   
    667             delta_i(ji,jj) = zdelta + rn_creepl 
     654            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 
     655            delta_i(ji,jj) = zdelta + rn_creepl * rswitch 
    668656 
    669657         END DO 
     
    678666 
    679667      !------------------------------------------------------------------------------! 
    680       ! 6) Control prints of residual and charge ellipse 
     668      ! 5) Control prints of residual and charge ellipse 
    681669      !------------------------------------------------------------------------------! 
    682670      ! 
     
    713701      ENDIF 
    714702      ! 
     703 
     704!!      ! velocity 
     705!!      IF ( iom_use( "uice_ipa" ) .OR. iom_use( "vice_ipa" ) .OR. iom_use( "icevel" ) ) THEN  
     706!!         CALL iom_put( "uice_ipa"     , u_ice      )       ! ice velocity u component 
     707!!         CALL iom_put( "vice_ipa"     , v_ice      )       ! ice velocity v component 
     708!!      ENDIF 
     709!!      IF ( iom_use( "iceconc_cat"  ) )  CALL iom_put( "iceconc_cat"      , a_i    )        ! area for categories 
     710!!      IF ( iom_use( "icethic_cat"  ) )  CALL iom_put( "icethic_cat"      , ht_i   )        ! thickness for categories 
     711      
     712       
    715713      CALL wrk_dealloc( jpi,jpj, z1_e1t0, z1_e2t0, zp_delt ) 
    716714      CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia ) 
    717715      CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV ) 
    718716      CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
    719       CALL wrk_dealloc( jpi,jpj, zswitch1, zswitch2, zfmask, zwf ) 
     717      CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf ) 
    720718 
    721719   END SUBROUTINE lim_rhg 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r6746 r6853  
    636636      INTEGER  :: ji, jk, jl             ! dummy loop indices 
    637637      INTEGER  :: ijpij, i_fill, jl0   
    638       REAL(wp) :: zarg, zV, zconv, zdh 
     638      REAL(wp) :: zarg, zV, zconv, zdh, zdv 
    639639      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables 
    640640      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables 
    641641      INTEGER , POINTER, DIMENSION(:)         ::   itest 
    642642  
    643       CALL wrk_alloc( 4, itest ) 
     643      Call wrk_alloc( 4, itest ) 
    644644      !-------------------------------------------------------------------- 
    645645      ! initialisation of variables 
     
    657657         IF( zhti(ji) > 0._wp ) THEN 
    658658 
    659          ! initialisation of tests 
    660          itest(:)  = 0 
     659            ! find which category (jl0) the input ice thickness falls into 
     660            jl0 = jpl 
     661            DO jl = 1, jpl 
     662               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 
     663                  jl0 = jl 
     664                  CYCLE 
     665               ENDIF 
     666            END DO 
     667 
     668            ! initialisation of tests 
     669            itest(:)  = 0 
    661670          
    662          i_fill = jpl + 1                                             !==================================== 
    663          DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories   
    664             ! iteration                                               !==================================== 
    665             i_fill = i_fill - 1 
     671            i_fill = jpl + 1                                             !==================================== 
     672            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories 
     673               ! iteration                                               !==================================== 
     674               i_fill = i_fill - 1 
     675                
     676               ! initialisation of ice variables for each try 
     677               zht_i(ji,1:jpl) = 0._wp 
     678               za_i (ji,1:jpl) = 0._wp 
     679               itest(:)        = 0            
     680                
     681               ! *** case very thin ice: fill only category 1 
     682               IF ( i_fill == 1 ) THEN 
     683                  zht_i(ji,1) = zhti(ji) 
     684                  za_i (ji,1) = zai (ji) 
     685                   
     686               ! *** case ice is thicker: fill categories >1 
     687               ELSE 
     688 
     689                  ! Fill ice thicknesses in the (i_fill-1) cat by hmean  
     690                  DO jl = 1, i_fill - 1 
     691                     zht_i(ji,jl) = hi_mean(jl) 
     692                  END DO 
     693                   
     694                  ! Concentrations in the (i_fill-1) categories  
     695                  za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl)) 
     696                  DO jl = 1, i_fill - 1 
     697                     IF ( jl /= jl0 ) THEN 
     698                        zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 
     699                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2) 
     700                     ENDIF 
     701                  END DO 
     702                   
     703                  ! Concentration in the last (i_fill) category 
     704                  za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) ) 
     705                   
     706                  ! Ice thickness in the last (i_fill) category 
     707                  zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) 
     708                  zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )  
     709                   
     710                  ! clem: correction if concentration of upper cat is greater than lower cat 
     711                  !       (it should be a gaussian around jl0 but sometimes it is not) 
     712                  IF ( jl0 /= jpl ) THEN 
     713                     DO jl = jpl, jl0+1, -1 
     714                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN 
     715                           zdv = zht_i(ji,jl) * za_i(ji,jl) 
     716                           zht_i(ji,jl    ) = 0._wp 
     717                           za_i (ji,jl    ) = 0._wp 
     718                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 ) 
     719                        END IF 
     720                     ENDDO 
     721                  ENDIF 
     722                
     723               ENDIF ! case ice is thick or thin 
     724                
     725               !--------------------- 
     726               ! Compatibility tests 
     727               !---------------------  
     728               ! Test 1: area conservation 
     729               zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) 
     730               IF ( zconv < epsi06 ) itest(1) = 1 
    666731             
    667             ! initialisation of ice variables for each try 
    668             zht_i(ji,1:jpl) = 0._wp 
    669             za_i (ji,1:jpl) = 0._wp 
    670              
    671             ! *** case very thin ice: fill only category 1 
    672             IF ( i_fill == 1 ) THEN 
    673                zht_i(ji,1) = zhti(ji) 
    674                za_i (ji,1) = zai (ji) 
    675  
    676             ! *** case ice is thicker: fill categories >1 
    677             ELSE 
    678  
    679                ! Fill ice thicknesses except the last one (i_fill) by hmean  
    680                DO jl = 1, i_fill - 1 
    681                   zht_i(ji,jl) = hi_mean(jl) 
    682                END DO 
     732               ! Test 2: volume conservation 
     733               zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) 
     734               IF ( zconv < epsi06 ) itest(2) = 1 
    683735                
    684                ! find which category (jl0) the input ice thickness falls into 
    685                jl0 = i_fill 
     736               ! Test 3: thickness of the last category is in-bounds ? 
     737               IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 
     738                
     739               ! Test 4: positivity of ice concentrations 
     740               itest(4) = 1 
    686741               DO jl = 1, i_fill 
    687                   IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 
    688                      jl0 = jl 
    689            CYCLE 
    690                   ENDIF 
    691                END DO 
    692                 
    693                ! Concentrations in the (i_fill-1) categories  
    694                za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl)) 
    695                DO jl = 1, i_fill - 1 
    696                   IF ( jl == jl0 ) CYCLE 
    697                   zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 
    698                   za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2) 
    699                END DO 
    700                 
    701                ! Concentration in the last (i_fill) category 
    702                za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) ) 
    703                 
    704                ! Ice thickness in the last (i_fill) category 
    705                zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) 
    706                zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill)  
    707                 
    708             ENDIF ! case ice is thick or thin 
    709              
    710             !--------------------- 
    711             ! Compatibility tests 
    712             !---------------------  
    713             ! Test 1: area conservation 
    714             zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) 
    715             IF ( zconv < epsi06 ) itest(1) = 1 
    716              
    717             ! Test 2: volume conservation 
    718             zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) 
    719             IF ( zconv < epsi06 ) itest(2) = 1 
    720              
    721             ! Test 3: thickness of the last category is in-bounds ? 
    722             IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 
    723              
    724             ! Test 4: positivity of ice concentrations 
    725             itest(4) = 1 
    726             DO jl = 1, i_fill 
    727                IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 
    728             END DO             
    729                                                            !============================ 
    730          END DO                                            ! end iteration on categories 
    731                                                            !============================ 
     742                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 
     743               END DO 
     744               !                                         !============================ 
     745            END DO                                       ! end iteration on categories 
     746               !                                         !============================ 
    732747         ENDIF ! if zhti > 0 
    733748      END DO ! i loop 
    734  
     749       
    735750      ! ------------------------------------------------ 
    736751      ! Adding Snow in each category where za_i is not 0 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r6764 r6853  
    9696         END DO 
    9797         CALL lbc_lnk( z2d, 'T', 1. ) 
    98          CALL iom_put( "uice_ipa"     , u_ice * zswi     )       ! ice velocity u component 
    99          CALL iom_put( "vice_ipa"     , v_ice * zswi     )       ! ice velocity v component 
    100          CALL iom_put( "icevel"       , z2d   * zswi     )       ! ice velocity module 
     98         CALL iom_put( "uice_ipa"     , u_ice      )       ! ice velocity u component 
     99         CALL iom_put( "vice_ipa"     , v_ice      )       ! ice velocity v component 
     100         CALL iom_put( "icevel"       , z2d        )       ! ice velocity module 
    101101      ENDIF 
    102102 
    103       IF ( iom_use( "tau_icebfr" ) )    CALL iom_put( "tau_icebfr"  , tau_icebfr * zswi      )  ! ice friction with ocean bottom (landfast ice)   
     103      IF ( iom_use( "tau_icebfr" ) )    CALL iom_put( "tau_icebfr"  , tau_icebfr             )  ! ice friction with ocean bottom (landfast ice)   
    104104      ! 
    105105      IF ( iom_use( "miceage" ) )       CALL iom_put( "miceage"     , om_i * zswi * z1_365   )  ! mean ice age 
     
    117117      CALL iom_put( "isnowhc"     , et_s  * zswi        )        ! snow total heat content 
    118118      CALL iom_put( "ibrinv"      , bvm_i * zswi * 100. )        ! brine volume 
    119       CALL iom_put( "utau_ice"    , utau_ice  * zswi    )        ! wind stress over ice along i-axis at I-point 
    120       CALL iom_put( "vtau_ice"    , vtau_ice  * zswi    )        ! wind stress over ice along j-axis at I-point 
     119      CALL iom_put( "utau_ice"    , utau_ice            )        ! wind stress over ice along i-axis at I-point 
     120      CALL iom_put( "vtau_ice"    , vtau_ice            )        ! wind stress over ice along j-axis at I-point 
    121121      CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation  
    122122      CALL iom_put( "micesalt"    , smt_i   * zswi      )        ! mean ice salinity 
    123123 
    124124      CALL iom_put( "icestr"      , strength * zswi )    ! ice strength 
    125       CALL iom_put( "idive"       , divu_i * 1.0e8   * zswi )    ! divergence 
    126       CALL iom_put( "ishear"      , shear_i * 1.0e8  * zswi )    ! shear 
     125      CALL iom_put( "idive"       , divu_i * 1.0e8     )    ! divergence 
     126      CALL iom_put( "ishear"      , shear_i * 1.0e8    )    ! shear 
    127127      CALL iom_put( "snowvol"     , vt_s   * zswi       )        ! snow volume 
    128128       
     
    232232      !!   4.0  !  2013-06  (C. Rousset) 
    233233      !!---------------------------------------------------------------------- 
    234       INTEGER, INTENT( in ) ::   kt               ! ocean time-step index) 
    235       INTEGER, INTENT( in ) ::   kid , kh_i        
     234      INTEGER, INTENT( in )   ::   kt               ! ocean time-step index) 
     235      INTEGER, INTENT( in )   ::   kid , kh_i 
     236      INTEGER                 ::   nz_i, jl 
     237      REAL(wp), DIMENSION(jpl) :: jcat 
    236238      !!---------------------------------------------------------------------- 
    237  
    238       CALL histdef( kid, "iicethic", "Ice thickness"           , "m"      ,   & 
    239       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    240       CALL histdef( kid, "iiceconc", "Ice concentration"       , "%"      ,   & 
    241       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    242       CALL histdef( kid, "iicetemp", "Ice temperature"         , "C"      ,   & 
    243       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    244       CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)"   , "m/s"    ,   & 
    245       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    246       CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)"   , "m/s"    ,   & 
    247       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    248       CALL histdef( kid, "iicestru", "i-Wind stress over ice (I-pt)", "Pa",   & 
    249       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    250       CALL histdef( kid, "iicestrv", "j-Wind stress over ice (I-pt)", "Pa",   & 
    251       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    252       CALL histdef( kid, "iicesflx", "Solar flux over ocean"     , "w/m2" ,   & 
    253       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    254       CALL histdef( kid, "iicenflx", "Non-solar flux over ocean" , "w/m2" ,   & 
     239      DO jl = 1, jpl 
     240         jcat(jl) = REAL(jl) 
     241      ENDDO 
     242       
     243      CALL histvert( kid, "ncatice", "Ice Categories","", jpl, jcat, nz_i, "up") 
     244 
     245      CALL histdef( kid, "sithic", "Ice thickness"           , "m"      ,   & 
     246      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     247      CALL histdef( kid, "siconc", "Ice concentration"       , "%"      ,   & 
     248      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     249      CALL histdef( kid, "sitemp", "Ice temperature"         , "C"      ,   & 
     250      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     251      CALL histdef( kid, "sivelu", "i-Ice speed "            , "m/s"    ,   & 
     252      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     253      CALL histdef( kid, "sivelv", "j-Ice speed "            , "m/s"    ,   & 
     254      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     255      CALL histdef( kid, "sistru", "i-Wind stress over ice " , "Pa"     ,   & 
     256      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     257      CALL histdef( kid, "sistrv", "j-Wind stress over ice " , "Pa"     ,   & 
     258      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     259      CALL histdef( kid, "sisflx", "Solar flux over ocean"     , "w/m2" ,   & 
     260      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     261      CALL histdef( kid, "sinflx", "Non-solar flux over ocean" , "w/m2" ,   & 
    255262      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    256263      CALL histdef( kid, "isnowpre", "Snow precipitation"      , "kg/m2/s",   & 
    257264      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    258       CALL histdef( kid, "iicesali", "Ice salinity"            , "PSU"    ,   & 
    259       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    260       CALL histdef( kid, "iicevolu", "Ice volume"              , "m"      ,   & 
    261       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    262       CALL histdef( kid, "iicedive", "Ice divergence"          , "10-8s-1",   & 
    263       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    264       CALL histdef( kid, "iicebopr", "Ice bottom production"   , "m/s"    ,   & 
    265       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    266       CALL histdef( kid, "iicedypr", "Ice dynamic production"  , "m/s"    ,   & 
    267       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    268       CALL histdef( kid, "iicelapr", "Ice open water prod"     , "m/s"    ,   & 
     265      CALL histdef( kid, "sisali", "Ice salinity"            , "PSU"    ,   & 
     266      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     267      CALL histdef( kid, "sivolu", "Ice volume"              , "m"      ,   & 
     268      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     269      CALL histdef( kid, "sidive", "Ice divergence"          , "10-8s-1",   & 
     270      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     271 
     272      CALL histdef( kid, "vfxbog", "Ice bottom production"   , "m/s"    ,   & 
     273      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     274      CALL histdef( kid, "vfxdyn", "Ice dynamic production"  , "m/s"    ,   & 
     275      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     276      CALL histdef( kid, "vfxopw", "Ice open water prod"     , "m/s"    ,   & 
    269277      &       jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    270       CALL histdef( kid, "iicesipr", "Snow ice production "    , "m/s"    ,   & 
    271       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    272       CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s"    ,   & 
    273       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    274       CALL histdef( kid, "iicebome", "Ice bottom melt"         , "m/s"    ,   & 
    275       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    276       CALL histdef( kid, "iicesume", "Ice surface melt"        , "m/s"    ,   & 
    277       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    278       CALL histdef( kid, "iisfxdyn", "Salt flux from dynmics"  , ""       ,   & 
    279       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    280       CALL histdef( kid, "iisfxres", "Salt flux from limupdate", ""       ,   & 
    281       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     278      CALL histdef( kid, "vfxsni", "Snow ice production "    , "m/s"    ,   & 
     279      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     280      CALL histdef( kid, "vfxres", "Ice prod from limupdate" , "m/s"    ,   & 
     281      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     282      CALL histdef( kid, "vfxbom", "Ice bottom melt"         , "m/s"    ,   & 
     283      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     284      CALL histdef( kid, "vfxsum", "Ice surface melt"        , "m/s"    ,   & 
     285      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     286 
     287      CALL histdef( kid, "sithicat", "Ice thickness"         , "m"      ,   & 
     288      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     289      CALL histdef( kid, "siconcat", "Ice concentration"     , "%"      ,   & 
     290      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     291      CALL histdef( kid, "sisalcat", "Ice salinity"           , ""      ,   & 
     292      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     293      CALL histdef( kid, "sitemcat", "Ice temperature"       , "C"      ,   & 
     294      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     295      CALL histdef( kid, "snthicat", "Snw thickness"         , "m"      ,   & 
     296      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     297      CALL histdef( kid, "sntemcat", "Snw temperature"       , "C"      ,   & 
     298      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
    282299 
    283300      CALL histend( kid, snc4set )   ! end of the file definition 
    284301 
    285       CALL histwrite( kid, "iicethic", kt, htm_i         , jpi*jpj, (/1/) )     
    286       CALL histwrite( kid, "iiceconc", kt, at_i          , jpi*jpj, (/1/) ) 
    287       CALL histwrite( kid, "iicetemp", kt, tm_i - rt0    , jpi*jpj, (/1/) ) 
    288       CALL histwrite( kid, "iicevelu", kt, u_ice          , jpi*jpj, (/1/) ) 
    289       CALL histwrite( kid, "iicevelv", kt, v_ice          , jpi*jpj, (/1/) ) 
    290       CALL histwrite( kid, "iicestru", kt, utau_ice       , jpi*jpj, (/1/) ) 
    291       CALL histwrite( kid, "iicestrv", kt, vtau_ice       , jpi*jpj, (/1/) ) 
    292       CALL histwrite( kid, "iicesflx", kt, qsr , jpi*jpj, (/1/) ) 
    293       CALL histwrite( kid, "iicenflx", kt, qns , jpi*jpj, (/1/) ) 
     302      CALL histwrite( kid, "sithic", kt, htm_i         , jpi*jpj, (/1/) )     
     303      CALL histwrite( kid, "siconc", kt, at_i          , jpi*jpj, (/1/) ) 
     304      CALL histwrite( kid, "sitemp", kt, tm_i - rt0    , jpi*jpj, (/1/) ) 
     305      CALL histwrite( kid, "sivelu", kt, u_ice          , jpi*jpj, (/1/) ) 
     306      CALL histwrite( kid, "sivelv", kt, v_ice          , jpi*jpj, (/1/) ) 
     307      CALL histwrite( kid, "sistru", kt, utau_ice       , jpi*jpj, (/1/) ) 
     308      CALL histwrite( kid, "sistrv", kt, vtau_ice       , jpi*jpj, (/1/) ) 
     309      CALL histwrite( kid, "sisflx", kt, qsr , jpi*jpj, (/1/) ) 
     310      CALL histwrite( kid, "sinflx", kt, qns , jpi*jpj, (/1/) ) 
    294311      CALL histwrite( kid, "isnowpre", kt, sprecip        , jpi*jpj, (/1/) ) 
    295       CALL histwrite( kid, "iicesali", kt, smt_i          , jpi*jpj, (/1/) ) 
    296       CALL histwrite( kid, "iicevolu", kt, vt_i           , jpi*jpj, (/1/) ) 
    297       CALL histwrite( kid, "iicedive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) ) 
    298  
    299       CALL histwrite( kid, "iicebopr", kt, wfx_bog        , jpi*jpj, (/1/) ) 
    300       CALL histwrite( kid, "iicedypr", kt, wfx_dyn        , jpi*jpj, (/1/) ) 
    301       CALL histwrite( kid, "iicelapr", kt, wfx_opw        , jpi*jpj, (/1/) ) 
    302       CALL histwrite( kid, "iicesipr", kt, wfx_sni        , jpi*jpj, (/1/) ) 
    303       CALL histwrite( kid, "iicerepr", kt, wfx_res        , jpi*jpj, (/1/) ) 
    304       CALL histwrite( kid, "iicebome", kt, wfx_bom        , jpi*jpj, (/1/) ) 
    305       CALL histwrite( kid, "iicesume", kt, wfx_sum        , jpi*jpj, (/1/) ) 
    306       CALL histwrite( kid, "iisfxdyn", kt, sfx_dyn        , jpi*jpj, (/1/) ) 
    307       CALL histwrite( kid, "iisfxres", kt, sfx_res        , jpi*jpj, (/1/) ) 
     312      CALL histwrite( kid, "sisali", kt, smt_i          , jpi*jpj, (/1/) ) 
     313      CALL histwrite( kid, "sivolu", kt, vt_i           , jpi*jpj, (/1/) ) 
     314      CALL histwrite( kid, "sidive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) ) 
     315 
     316      CALL histwrite( kid, "vfxbog", kt, wfx_bog        , jpi*jpj, (/1/) ) 
     317      CALL histwrite( kid, "vfxdyn", kt, wfx_dyn        , jpi*jpj, (/1/) ) 
     318      CALL histwrite( kid, "vfxopw", kt, wfx_opw        , jpi*jpj, (/1/) ) 
     319      CALL histwrite( kid, "vfxsni", kt, wfx_sni        , jpi*jpj, (/1/) ) 
     320      CALL histwrite( kid, "vfxres", kt, wfx_res        , jpi*jpj, (/1/) ) 
     321      CALL histwrite( kid, "vfxbom", kt, wfx_bom        , jpi*jpj, (/1/) ) 
     322      CALL histwrite( kid, "vfxsum", kt, wfx_sum        , jpi*jpj, (/1/) ) 
     323 
     324      CALL histwrite( kid, "sithicat", kt, ht_i        , jpi*jpj*jpl, (/1/) )     
     325      CALL histwrite( kid, "siconcat", kt, a_i         , jpi*jpj*jpl, (/1/) )     
     326      CALL histwrite( kid, "sisalcat", kt, sm_i        , jpi*jpj*jpl, (/1/) )     
     327      CALL histwrite( kid, "sitemcat", kt, tm_i - rt0  , jpi*jpj*jpl, (/1/) )     
     328      CALL histwrite( kid, "snthicat", kt, ht_s        , jpi*jpj*jpl, (/1/) )     
     329      CALL histwrite( kid, "sntemcat", kt, tm_su - rt0 , jpi*jpj*jpl, (/1/) )     
    308330 
    309331      ! Close the file 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90

    r4205 r6853  
    4343   INTEGER, PARAMETER, PUBLIC ::   jp_i1    = 204      !: write INTEGER(1) 
    4444 
    45    INTEGER, PARAMETER, PUBLIC ::   jpmax_files  = 100   !: maximum number of simultaneously opened file 
    46    INTEGER, PARAMETER, PUBLIC ::   jpmax_vars   = 600 !: maximum number of variables in one file 
     45   INTEGER, PARAMETER, PUBLIC ::   jpmax_files  = 100  !: maximum number of simultaneously opened file 
     46   INTEGER, PARAMETER, PUBLIC ::   jpmax_vars   = 1200 !: maximum number of variables in one file 
    4747   INTEGER, PARAMETER, PUBLIC ::   jpmax_dims   =  4   !: maximum number of dimensions for one variable 
    4848   INTEGER, PARAMETER, PUBLIC ::   jpmax_digits =  5   !: maximum number of digits for the cpu number in the file name 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r6824 r6853  
    947947 
    948948 
     949#if defined key_lim3 
    949950   SUBROUTINE Cdn10_Lupkes2012( Cd ) 
    950951      !!---------------------------------------------------------------------- 
     
    991992       
    992993   END SUBROUTINE Cdn10_Lupkes2012 
    993        
     994#endif 
     995    
    994996   !!====================================================================== 
    995997END MODULE sbcblk_core 
Note: See TracChangeset for help on using the changeset viewer.