Ignore:
Timestamp:
2017-09-12T20:46:13+02:00 (3 years ago)
Author:
clem
Message:

changes in style - part6 - one more round

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90

    r8516 r8517  
    3737   PUBLIC   ice_adv_init   ! called by icestp 
    3838 
     39   INTEGER ::              nice_adv   ! choice of the type of advection scheme 
     40   !                                        ! associated indices: 
     41   INTEGER, PARAMETER ::   np_advPRA = 1   ! Prather scheme 
     42   INTEGER, PARAMETER ::   np_advUMx = 2   ! Ultimate-Macho scheme 
     43 
    3944   !! * Substitution 
    4045#  include "vectopt_loop_substitute.h90" 
     
    6166      !!---------------------------------------------------------------------- 
    6267      INTEGER, INTENT(in) ::   kt   ! number of iteration 
    63       ! 
    64       INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices 
    65       INTEGER  ::   initad                  ! number of sub-timestep for the advection 
    66       REAL(wp) ::   zcfl , zusnit           !   -      - 
    67       CHARACTER(len=80) :: cltmp 
    68       ! 
    69       REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
    70       REAL(wp) ::    zdv 
    71       REAL(wp), DIMENSION(jpi,jpj)           ::   zatold, zeiold, zesold, zsmvold  
    72       REAL(wp), DIMENSION(jpi,jpj,jpl)       ::   zhimax, zviold, zvsold 
    7368      !!--------------------------------------------------------------------- 
    7469      IF( nn_timing == 1 )  CALL timing_start('iceadv') 
     
    8075      ENDIF 
    8176       
    82       CALL ice_var_agg( 1 ) ! integrated values + ato_i 
    83  
    8477      ! conservation test 
    85       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    86        
    87       ! store old values for diag 
    88       zviold (:,:,:) = v_i(:,:,:) 
    89       zvsold (:,:,:) = v_s(:,:,:) 
    90       zsmvold(:,:)   = SUM( smv_i(:,:,:), dim=3 ) 
    91       zeiold (:,:)   = et_i(:,:) 
    92       zesold (:,:)   = et_s(:,:) 
    93  
    94       ! Thickness correction init. 
    95       zatold(:,:) = at_i 
    96       WHERE( a_i(:,:,:) >= epsi20 ) 
    97          ht_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:) 
    98          ht_s(:,:,:) = v_s(:,:,:) / a_i(:,:,:) 
    99       ELSEWHERE 
    100          ht_i(:,:,:) = 0._wp 
    101          ht_s(:,:,:) = 0._wp          
    102       END WHERE 
    103           
    104       ! Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick 
    105       zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:) 
    106       DO jl = 1, jpl 
    107          DO jj = 2, jpjm1 
    108             DO ji = 2, jpim1 
    109 !!gm use of MAXVAL here is very probably less efficient than expending the 9 values 
    110                zhimax(ji,jj,jl) = MAX( epsi20, MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) ) ) 
    111             END DO 
    112          END DO 
    113       END DO 
    114       CALL lbc_lnk( zhimax(:,:,:), 'T', 1. ) 
    115        
     78      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceadv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     79                      
    11680      !---------- 
    11781      ! Advection 
    11882      !---------- 
    119       IF( ln_adv_UMx ) THEN                    !-- ULTIMATE-MACHO scheme 
     83      SELECT CASE( nice_adv ) 
     84      !                                !-----------------------! 
     85      CASE( np_advUMx )                ! ULTIMATE-MACHO scheme ! 
     86         !                             !-----------------------! 
    12087         CALL ice_adv_umx( kt, u_ice, v_ice,  & 
    12188            &              ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    122           
    123       ELSEIF( ln_adv_Pra ) THEN                !-- PRATHER scheme 
     89      !                                !-----------------------! 
     90      CASE( np_advPRA )                ! PRATHER scheme        ! 
     91         !                             !-----------------------! 
    12492         CALL ice_adv_prather( kt, u_ice, v_ice,  & 
    12593            &                  ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    126           
    127       ENDIF 
     94      END SELECT 
    12895 
    129       ! total ice fraction 
    130       at_i(:,:) = a_i(:,:,1) 
    131       DO jl = 2, jpl 
    132          at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    133       END DO 
    134        
    13596      !------------ 
    13697      ! diagnostics 
    13798      !------------ 
    138       DO jj = 1, jpj 
    139          DO ji = 1, jpi 
    140             diag_trp_ei (ji,jj) = ( SUM( e_i  (ji,jj,1:nlay_i,:) ) -  zeiold(ji,jj) ) * r1_rdtice 
    141             diag_trp_es (ji,jj) = ( SUM( e_s  (ji,jj,1:nlay_s,:) ) -  zesold(ji,jj) ) * r1_rdtice 
    142             diag_trp_smv(ji,jj) = ( SUM( smv_i(ji,jj,:)          ) - zsmvold(ji,jj) ) * r1_rdtice 
    143             diag_trp_vi (ji,jj) =   SUM(   v_i(ji,jj,:)            -  zviold(ji,jj,:) ) * r1_rdtice 
    144             diag_trp_vs (ji,jj) =   SUM(   v_s(ji,jj,:)            -  zvsold(ji,jj,:) ) * r1_rdtice 
    145          END DO 
    146       END DO 
     99      diag_trp_ei (:,:) = SUM(SUM( e_i  (:,:,1:nlay_i,:) - e_i_b  (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_rdtice 
     100      diag_trp_es (:,:) = SUM(SUM( e_s  (:,:,1:nlay_s,:) - e_s_b  (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_rdtice 
     101      diag_trp_smv(:,:) = SUM(     smv_i(:,:,:)          - smv_i_b(:,:,:)                  , dim=3 ) * r1_rdtice 
     102      diag_trp_vi (:,:) = SUM(     v_i  (:,:,:)          - v_i_b  (:,:,:)                  , dim=3 ) * r1_rdtice 
     103      diag_trp_vs (:,:) = SUM(     v_s  (:,:,:)          - v_s_b  (:,:,:)                  , dim=3 ) * r1_rdtice 
    147104      IF( iom_use('icetrp') )   CALL iom_put( "icetrp" , diag_trp_vi * rday  )         ! ice volume transport 
    148105      IF( iom_use('snwtrp') )   CALL iom_put( "snwtrp" , diag_trp_vs * rday  )         ! snw volume transport 
     
    150107      IF( iom_use('deitrp') )   CALL iom_put( "deitrp" , diag_trp_ei         )         ! advected ice enthalpy (W/m2) 
    151108      IF( iom_use('destrp') )   CALL iom_put( "destrp" , diag_trp_es         )         ! advected snw enthalpy (W/m2) 
    152        
    153       !-------------------------------------- 
    154       ! Thickness correction in case too high 
    155       !-------------------------------------- 
    156       IF( nn_icedyn == 2 ) THEN 
    157          ! 
    158          CALL ice_var_zapsmall                       !-- zap small areas 
    159          ! 
    160          DO jl = 1, jpl 
    161             DO jj = 1, jpj 
    162                DO ji = 1, jpi 
    163                   IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to zhimax 
    164                      ! 
    165                      ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 
    166                      ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / a_i(ji,jj,jl) 
    167                      zdv = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl)   
    168                      ! 
    169                      IF ( ( zdv >  0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. & 
    170                         & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN 
    171                         a_i (ji,jj,jl) = ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / zhimax(ji,jj,jl) 
    172                         ht_i(ji,jj,jl) =   v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    173                      ENDIF 
    174                      ! 
    175                   ENDIF 
    176                END DO 
    177             END DO 
    178          END DO 
    179                    
    180          WHERE( ht_i(:,:,jpl) > hi_max(jpl) )        !-- bound ht_i to hi_max (99 m) 
    181             ht_i(:,:,jpl) = hi_max(jpl) 
    182             a_i (:,:,jpl) = v_i(:,:,jpl) / hi_max(jpl) 
    183          END WHERE 
    184  
    185          IF ( nn_pnd_scheme > 0 ) THEN               !-- correct pond fraction to avoid a_ip > a_i 
    186             WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:) 
    187          ENDIF 
    188          ! 
    189       ENDIF 
    190           
    191       !------------------------------------------------------------ 
    192       ! Impose a_i < amax if no ridging/rafting or in mono-category 
    193       !------------------------------------------------------------ 
    194       IF( l_piling ) THEN                            !-- simple conservative piling, comparable with 1-cat models 
    195          at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    196          DO jl = 1, jpl 
    197             WHERE( at_i(:,:) > epsi20 ) 
    198                a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  ) 
    199             END WHERE 
    200          END DO 
    201       ENDIF 
    202        
    203       ! agglomerate variables  
    204       vt_i(:,:) = SUM( v_i(:,:,:), dim=3 ) 
    205       vt_s(:,:) = SUM( v_s(:,:,:), dim=3 ) 
    206       at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    207  
    208       ! MV MP 2016 (remove once we get rid of a_i_frac and ht_i) 
    209       IF ( nn_pnd_scheme > 0 ) THEN 
    210           at_ip(:,:) = SUM( a_ip(:,:,:), dim = 3 ) 
    211           vt_ip(:,:) = SUM( v_ip(:,:,:), dim = 3 ) 
    212       ENDIF 
    213       ! END MP 2016 
    214        
    215       ! open water = 1 if at_i=0  
    216       WHERE( at_i == 0._wp ) ato_i = 1._wp  
    217        
    218       ! conservation test 
    219       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    220          
    221       ! -------------- 
    222       ! control prints 
    223       ! -------------- 
    224       IF( ln_icectl )   CALL ice_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' ) 
    225       ! 
    226       IF( nn_timing == 1 )  CALL timing_stop('iceadv') 
     109                
     110      ! controls 
     111      IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'iceadv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     112      IF( ln_icectl      )   CALL ice_prt     (kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ')                       ! prints 
     113      IF( nn_timing == 1 )   CALL timing_stop ('iceadv')                                                             ! timing 
    227114      ! 
    228115   END SUBROUTINE ice_adv 
     
    241128      !! ** input   :   Namelist namice_adv 
    242129      !!------------------------------------------------------------------- 
    243       INTEGER ::   ios   ! Local integer output status for namelist read 
     130      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
    244131      !! 
    245132      NAMELIST/namice_adv/ ln_adv_Pra, ln_adv_UMx, nn_UMx 
     
    260147         WRITE(numout,*) '~~~~~~~~~~~~' 
    261148         WRITE(numout,*) '   Namelist namice_adv' 
    262          WRITE(numout,*) '      advection scheme for ice transport (limtrp)' 
    263          WRITE(numout,*) '         type of advection scheme (Prather)                     ln_adv_Pra = ', ln_adv_Pra  
    264          WRITE(numout,*) '         type of advection scheme (Ulimate-Macho)               ln_adv_UMx = ', ln_adv_UMx  
    265          WRITE(numout,*) '            order of the Ultimate-Macho scheme                      nn_UMx = ', nn_UMx 
     149         WRITE(numout,*) '      type of advection scheme (Prather)                     ln_adv_Pra = ', ln_adv_Pra  
     150         WRITE(numout,*) '      type of advection scheme (Ulimate-Macho)               ln_adv_UMx = ', ln_adv_UMx  
     151         WRITE(numout,*) '         order of the Ultimate-Macho scheme                      nn_UMx = ', nn_UMx 
    266152      ENDIF 
    267153      ! 
    268       IF ( ( ln_adv_Pra .AND. ln_adv_UMx ) .OR. ( .NOT.ln_adv_Pra .AND. .NOT.ln_adv_UMx ) ) THEN 
    269          CALL ctl_stop( 'ice_adv_init: choose one and only one ice advection scheme (ln_adv_Pra or ln_adv_UMx)' ) 
    270       ENDIF 
     154      !                             !== set the choice of ice advection ==! 
     155      ioptio = 0  
     156      IF( ln_adv_Pra ) THEN   ;   ioptio = ioptio + 1   ;   nice_adv = np_advPRA    ;   ENDIF 
     157      IF( ln_adv_UMx ) THEN   ;   ioptio = ioptio + 1   ;   nice_adv = np_advUMx    ;   ENDIF 
     158      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_adv_init: choose one and only one ice advection scheme (ln_adv_Pra or ln_adv_UMx)' ) 
    271159      ! 
    272160   END SUBROUTINE ice_adv_init 
Note: See TracChangeset for help on using the changeset viewer.