Ignore:
Timestamp:
2015-07-15T17:46:12+02:00 (5 years ago)
Author:
andrewryan
Message:

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90

    r5034 r5600  
    1313   !!---------------------------------------------------------------------- 
    1414   !!   lim_hdf       : diffusion trend on sea-ice variable 
     15   !!   lim_hdf_init  : initialisation of diffusion trend on sea-ice variable 
    1516   !!---------------------------------------------------------------------- 
    1617   USE dom_oce        ! ocean domain 
     
    2627   PRIVATE 
    2728 
    28    PUBLIC   lim_hdf     ! called by lim_tra 
    29  
    30    LOGICAL  ::   linit = .TRUE.              ! initialization flag (set to flase after the 1st call) 
    31    REAL(wp) ::   epsi04 = 1.e-04              ! constant 
     29   PUBLIC   lim_hdf         ! called by lim_trp 
     30   PUBLIC   lim_hdf_init    ! called by sbc_lim_init 
     31 
     32   LOGICAL  ::   linit = .TRUE.                             ! initialization flag (set to flase after the 1st call) 
     33   INTEGER  ::   nn_convfrq                                 !:  convergence check frequency of the Crant-Nicholson scheme 
    3234   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   efact   ! metric coefficient 
    3335 
     
    5456      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   ptab    ! Field on which the diffusion is applied 
    5557      ! 
    56       INTEGER  ::  ji, jj                   ! dummy loop indices 
    57       INTEGER  ::  its, iter, ierr          ! local integers 
    58       REAL(wp) ::   zalfa, zrlxint, zconv   ! local scalars 
    59       REAL(wp), POINTER, DIMENSION(:,:) ::   zrlx, zflu, zflv, zdiv0, zdiv, ztab0 
    60       CHARACTER(lc) ::   charout   ! local character 
     58      INTEGER                           ::  ji, jj                    ! dummy loop indices 
     59      INTEGER                           ::  iter, ierr           ! local integers 
     60      REAL(wp)                          ::  zrlxint, zconv     ! local scalars 
     61      REAL(wp), POINTER, DIMENSION(:,:) ::  zrlx, zflu, zflv, zdiv0, zdiv, ztab0 
     62      CHARACTER(lc)                     ::  charout                   ! local character 
     63      REAL(wp), PARAMETER               ::  zrelax = 0.5_wp           ! relaxation constant for iterative procedure 
     64      REAL(wp), PARAMETER               ::  zalfa  = 0.5_wp           ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit 
     65      INTEGER , PARAMETER               ::  its    = 100              ! Maximum number of iteration 
    6166      !!------------------------------------------------------------------- 
    6267       
     
    7176         DO jj = 2, jpjm1   
    7277            DO ji = fs_2 , fs_jpim1   ! vector opt. 
    73                efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) ) 
     78               efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) * r1_e12t(ji,jj) 
    7479            END DO 
    7580         END DO 
     
    7782      ENDIF 
    7883      !                             ! Time integration parameters 
    79       zalfa = 0.5_wp                      ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit 
    80       its   = 100                         ! Maximum number of iteration 
    8184      ! 
    8285      ztab0(:, : ) = ptab(:,:)      ! Arrays initialization 
     
    9194      iter  = 0 
    9295      ! 
    93       DO WHILE( zconv > ( 2._wp * epsi04 ) .AND. iter <= its )   ! Sub-time step loop 
     96      DO WHILE( zconv > ( 2._wp * 1.e-04 ) .AND. iter <= its )   ! Sub-time step loop 
    9497         ! 
    9598         iter = iter + 1                                 ! incrementation of the sub-time step number 
     
    97100         DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction 
    98101            DO ji = 1 , fs_jpim1   ! vector opt. 
    99                zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 
    100                zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) / e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 
     102               zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 
     103               zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 
    101104            END DO 
    102105         END DO 
     
    104107         DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes 
    105108            DO ji = fs_2 , fs_jpim1   ! vector opt.  
    106                zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   & 
    107                   &            + zflv(ji,jj) - zflv(ji  ,jj-1)  ) / ( e1t (ji,jj) * e2t (ji,jj) ) 
     109               zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj) 
    108110            END DO 
    109111         END DO 
     
    115117               zrlxint = (   ztab0(ji,jj)    & 
    116118                  &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj) )   & 
    117                   &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )  )                             &  
    118                   &    / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) ) 
    119                zrlx(ji,jj) = ptab(ji,jj) + om * ( zrlxint - ptab(ji,jj) ) 
     119                  &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )                               &  
     120                  &      ) / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) ) 
     121               zrlx(ji,jj) = ptab(ji,jj) + zrelax * ( zrlxint - ptab(ji,jj) ) 
    120122            END DO 
    121123         END DO 
    122124         CALL lbc_lnk( zrlx, 'T', 1. )                   ! lateral boundary condition 
    123125         ! 
    124          zconv = 0._wp                                   ! convergence test 
    125          DO jj = 2, jpjm1 
    126             DO ji = fs_2, fs_jpim1 
    127                zconv = MAX( zconv, ABS( zrlx(ji,jj) - ptab(ji,jj) )  ) 
    128             END DO 
    129          END DO 
    130          IF( lk_mpp )   CALL mpp_max( zconv )            ! max over the global domain 
     126         IF ( MOD( iter, nn_convfrq ) == 0 )  THEN    ! convergence test every nn_convfrq iterations (perf. optimization) 
     127            zconv = 0._wp 
     128            DO jj = 2, jpjm1 
     129               DO ji = fs_2, fs_jpim1 
     130                  zconv = MAX( zconv, ABS( zrlx(ji,jj) - ptab(ji,jj) )  ) 
     131               END DO 
     132            END DO 
     133            IF( lk_mpp )   CALL mpp_max( zconv )      ! max over the global domain 
     134         ENDIF 
    131135         ! 
    132136         ptab(:,:) = zrlx(:,:) 
     
    138142      DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction 
    139143         DO ji = 1 , fs_jpim1   ! vector opt. 
    140             zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 
    141             zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) / e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 
     144            zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 
     145            zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 
    142146         END DO 
    143147      END DO 
     
    145149      DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes 
    146150         DO ji = fs_2 , fs_jpim1   ! vector opt.  
    147             zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   & 
    148                  &            + zflv(ji,jj) - zflv(ji  ,jj-1)  ) / ( e1t (ji,jj) * e2t (ji,jj) ) 
     151            zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj) 
    149152            ptab(ji,jj) = ztab0(ji,jj) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj) ) 
    150153         END DO 
     
    164167   END SUBROUTINE lim_hdf 
    165168 
     169    
     170   SUBROUTINE lim_hdf_init 
     171      !!------------------------------------------------------------------- 
     172      !!                  ***  ROUTINE lim_hdf_init  *** 
     173      !! 
     174      !! ** Purpose : Initialisation of horizontal diffusion of sea-ice  
     175      !! 
     176      !! ** Method  : Read the namicehdf namelist 
     177      !! 
     178      !! ** input   : Namelist namicehdf 
     179      !!------------------------------------------------------------------- 
     180      INTEGER  ::   ios                 ! Local integer output status for namelist read 
     181      NAMELIST/namicehdf/ nn_convfrq 
     182      !!------------------------------------------------------------------- 
     183      ! 
     184      IF(lwp) THEN 
     185         WRITE(numout,*) 
     186         WRITE(numout,*) 'lim_hdf : Ice horizontal diffusion' 
     187         WRITE(numout,*) '~~~~~~~' 
     188      ENDIF 
     189      ! 
     190      REWIND( numnam_ice_ref )              ! Namelist namicehdf in reference namelist : Ice horizontal diffusion 
     191      READ  ( numnam_ice_ref, namicehdf, IOSTAT = ios, ERR = 901) 
     192901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicehdf in reference namelist', lwp ) 
     193 
     194      REWIND( numnam_ice_cfg )              ! Namelist namicehdf in configuration namelist : Ice horizontal diffusion 
     195      READ  ( numnam_ice_cfg, namicehdf, IOSTAT = ios, ERR = 902 ) 
     196902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicehdf in configuration namelist', lwp ) 
     197      IF(lwm) WRITE ( numoni, namicehdf ) 
     198      ! 
     199      IF(lwp) THEN                          ! control print 
     200         WRITE(numout,*) 
     201         WRITE(numout,*)'   Namelist of ice parameters for ice horizontal diffusion computation ' 
     202         WRITE(numout,*)'      convergence check frequency of the Crant-Nicholson scheme   nn_convfrq   = ', nn_convfrq 
     203      ENDIF 
     204      ! 
     205   END SUBROUTINE lim_hdf_init 
    166206#else 
    167207   !!---------------------------------------------------------------------- 
    168208   !!   Default option          Dummy module           NO LIM sea-ice model 
    169209   !!---------------------------------------------------------------------- 
    170 CONTAINS 
    171    SUBROUTINE lim_hdf         ! Empty routine 
    172    END SUBROUTINE lim_hdf 
    173210#endif 
    174211 
Note: See TracChangeset for help on using the changeset viewer.