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 6126 for branches/2015/dev_r5546_CNRS19_HPC_scalability/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90 – NEMO

Ignore:
Timestamp:
2015-12-18T13:58:27+01:00 (8 years ago)
Author:
mcastril
Message:

lim_hdf routine substituted by lim_hdf_multiple, all limtrp calls use the multiple version

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5546_CNRS19_HPC_scalability/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90

    r6052 r6126  
    77   !!             -   !  2001-05 (G. Madec, R. Hordoir) opa norm 
    88   !!            1.0  !  2002-08 (C. Ethe)  F90, free form 
    9    !!            3.0  !  2015-08 (O. Tintó and M. Castrillo)  added lim_hdf_multiple 
     9   !!            3.0  !  2015-08 (O. Tintó and M. Castrillo)  added lim_hdf (multiple) 
    1010   !!---------------------------------------------------------------------- 
    1111#if defined key_lim3 
     
    2828   PRIVATE 
    2929 
    30    PUBLIC   lim_hdf         ! called by lim_trp 
    31    PUBLIC   lim_hdf_multiple ! called by lim_trp 
     30   PUBLIC   lim_hdf ! called by lim_trp 
    3231   PUBLIC   lim_hdf_init    ! called by sbc_lim_init 
    3332 
     
    4544CONTAINS 
    4645 
    47    SUBROUTINE lim_hdf( ptab ) 
     46   SUBROUTINE lim_hdf( ptab , ihdf_vars , jpl , nlay_i ) 
    4847      !!------------------------------------------------------------------- 
    4948      !!                  ***  ROUTINE lim_hdf  *** 
     
    5655      !! ** Action  :    update ptab with the diffusive contribution 
    5756      !!------------------------------------------------------------------- 
    58       REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   ptab    ! Field on which the diffusion is applied 
    59       ! 
    60       INTEGER                           ::  ji, jj                    ! dummy loop indices 
    61       INTEGER                           ::  iter, ierr           ! local integers 
    62       REAL(wp)                          ::  zrlxint, zconv     ! local scalars 
    63       REAL(wp), POINTER, DIMENSION(:,:) ::  zrlx, zflu, zflv, zdiv0, zdiv, ztab0 
    64       CHARACTER(lc)                     ::  charout                   ! local character 
    65       REAL(wp), PARAMETER               ::  zrelax = 0.5_wp           ! relaxation constant for iterative procedure 
    66       REAL(wp), PARAMETER               ::  zalfa  = 0.5_wp           ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit 
    67       INTEGER , PARAMETER               ::  its    = 100              ! Maximum number of iteration 
    68       !!------------------------------------------------------------------- 
    69        
    70       CALL wrk_alloc( jpi, jpj, zrlx, zflu, zflv, zdiv0, zdiv, ztab0 ) 
    71  
    72       !                       !==  Initialisation  ==! 
    73       ! 
    74       IF( linit ) THEN              ! Metric coefficient (compute at the first call and saved in efact) 
    75          ALLOCATE( efact(jpi,jpj) , STAT=ierr ) 
    76          IF( lk_mpp    )   CALL mpp_sum( ierr ) 
    77          IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'lim_hdf : unable to allocate arrays' ) 
    78          DO jj = 2, jpjm1   
    79             DO ji = fs_2 , fs_jpim1   ! vector opt. 
    80                efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) * r1_e12t(ji,jj) 
    81             END DO 
    82          END DO 
    83          linit = .FALSE. 
    84       ENDIF 
    85       !                             ! Time integration parameters 
    86       ! 
    87       ztab0(:, : ) = ptab(:,:)      ! Arrays initialization 
    88       zdiv0(:, 1 ) = 0._wp 
    89       zdiv0(:,jpj) = 0._wp 
    90       zflu (jpi,:) = 0._wp    
    91       zflv (jpi,:) = 0._wp 
    92       zdiv0(1,  :) = 0._wp 
    93       zdiv0(jpi,:) = 0._wp 
    94  
    95       zconv = 1._wp           !==  horizontal diffusion using a Crant-Nicholson scheme  ==! 
    96       iter  = 0 
    97       ! 
    98       DO WHILE( zconv > ( 2._wp * 1.e-04 ) .AND. iter <= its )   ! Sub-time step loop 
    99          ! 
    100          iter = iter + 1                                 ! incrementation of the sub-time step number 
    101          ! 
    102          DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction 
    103             DO ji = 1 , fs_jpim1   ! vector opt. 
    104                zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 
    105                zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 
    106             END DO 
    107          END DO 
    108          ! 
    109          DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes 
    110             DO ji = fs_2 , fs_jpim1   ! vector opt.  
    111                zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj) 
    112             END DO 
    113          END DO 
    114          ! 
    115          IF( iter == 1 )   zdiv0(:,:) = zdiv(:,:)        ! save the 1st evaluation of the diffusive trend in zdiv0 
    116          ! 
    117          DO jj = 2, jpjm1                                ! iterative evaluation 
    118             DO ji = fs_2 , fs_jpim1   ! vector opt. 
    119                zrlxint = (   ztab0(ji,jj)    & 
    120                   &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj) )   & 
    121                   &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )                               &  
    122                   &      ) / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) ) 
    123                zrlx(ji,jj) = ptab(ji,jj) + zrelax * ( zrlxint - ptab(ji,jj) ) 
    124             END DO 
    125          END DO 
    126          CALL lbc_lnk( zrlx, 'T', 1. )                   ! lateral boundary condition 
    127          ! 
    128          IF ( MOD( iter - 1 , nn_convfrq ) == 0 )  THEN    ! convergence test every nn_convfrq iterations (perf. optimization) 
    129             zconv = 0._wp 
    130             DO jj = 2, jpjm1 
    131                DO ji = fs_2, fs_jpim1 
    132                   zconv = MAX( zconv, ABS( zrlx(ji,jj) - ptab(ji,jj) )  ) 
    133                END DO 
    134             END DO 
    135             IF( lk_mpp )   CALL mpp_max( zconv )      ! max over the global domain 
    136          ENDIF 
    137          ! 
    138          ptab(:,:) = zrlx(:,:) 
    139          ! 
    140       END DO                                       ! end of sub-time step loop 
    141  
    142       ! ----------------------- 
    143       !!! final step (clem) !!! 
    144       DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction 
    145          DO ji = 1 , fs_jpim1   ! vector opt. 
    146             zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 
    147             zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 
    148          END DO 
    149       END DO 
    150       ! 
    151       DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes 
    152          DO ji = fs_2 , fs_jpim1   ! vector opt.  
    153             zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj) 
    154             ptab(ji,jj) = ztab0(ji,jj) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj) ) 
    155          END DO 
    156       END DO 
    157       CALL lbc_lnk( ptab, 'T', 1. )                   ! lateral boundary condition 
    158       !!! final step (clem) !!! 
    159       ! ----------------------- 
    160  
    161       IF(ln_ctl)   THEN 
    162          zrlx(:,:) = ptab(:,:) - ztab0(:,:) 
    163          WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter 
    164          CALL prt_ctl( tab2d_1=zrlx, clinfo1=charout ) 
    165       ENDIF 
    166       ! 
    167       CALL wrk_dealloc( jpi, jpj, zrlx, zflu, zflv, zdiv0, zdiv, ztab0 ) 
    168       ! 
    169    END SUBROUTINE lim_hdf 
    170     
    171  
    172    SUBROUTINE lim_hdf_multiple( ptab , ihdf_vars , jpl , nlay_i ) 
    173       !!------------------------------------------------------------------- 
    174       !!                  ***  ROUTINE lim_hdf  *** 
    175       !! 
    176       !! ** purpose :   Compute and add the diffusive trend on sea-ice variables 
    177       !! 
    178       !! ** method  :   Second order diffusive operator evaluated using a 
    179       !!              Cranck-Nicholson time Scheme. 
    180       !! 
    181       !! ** Action  :    update ptab with the diffusive contribution 
    182       !!------------------------------------------------------------------- 
    18357      INTEGER                           :: jpl, nlay_i, isize, ihdf_vars 
    18458      REAL(wp),  DIMENSION(:,:,:), INTENT( inout ),TARGET ::   ptab    ! Field on which the diffusion is applied 
    185       REAL(wp), POINTER, DIMENSION(:,:,:)        ::   pahu3D , pahv3D 
    18659      ! 
    18760      INTEGER                           ::  ji, jj, jk, jl , jm               ! dummy loop indices 
     
    20477 
    20578      !                       !==  Initialisation  ==! 
    206          isize = jpl*(ihdf_vars+nlay_i) 
     79      ! +1 open water diffusion 
     80      isize = jpl*(ihdf_vars+nlay_i)+1 
    20781      ALLOCATE( zconv (isize) ) 
    20882      ALLOCATE( pt2d_array(isize) , zrlx_array(isize) ) 
     
    21286      CALL wrk_alloc( jpi, jpj, isize, zrlx, zdiv0, ztab0 ) 
    21387      CALL wrk_alloc( jpi, jpj, zflu, zflv, zdiv ) 
    214       CALL wrk_alloc( jpi, jpj, jpl, pahu3D , pahv3D ) 
    215  
    216  
    217       DO jl = 1 , jpl 
    218          jm = (jl-1)*(ihdf_vars+nlay_i)+1 
    219          DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row 
    220             DO ji = 1 , fs_jpim1   ! vector opt. 
    221                pahu3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ptab(ji  ,jj,jm) ) ) )   & 
    222                &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ptab(ji+1, jj, jm ) ) ) ) * ahiu(ji,jj) 
    223                pahv3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ptab(ji, jj, jm ) ) ) )   & 
    224                &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- ptab(ji, jj+1, jm ) ) ) ) * ahiv(ji,jj) 
    225             END DO 
    226          END DO 
    227       END DO 
    22888 
    22989      DO jk= 1 , isize 
     
    266126         iter = iter + 1                                 ! incrementation of the sub-time step number 
    267127         ! 
    268  
    269128         DO jk = 1 , isize 
    270129            jl = (jk-1) /( ihdf_vars+nlay_i)+1 
     
    353212      CALL wrk_dealloc( jpi, jpj, isize, zrlx, zdiv0, ztab0 ) 
    354213      CALL wrk_dealloc( jpi, jpj, zflu, zflv, zdiv ) 
    355       CALL wrk_dealloc( jpi, jpj, jpl, pahu3D , pahv3D ) 
    356214 
    357215      DEALLOCATE( zconv ) 
     
    360218      DEALLOCATE( psgn_array ) 
    361219      ! 
    362    END SUBROUTINE lim_hdf_multiple 
     220   END SUBROUTINE lim_hdf 
    363221 
    364222 
     
    408266   !!====================================================================== 
    409267END MODULE limhdf 
     268 
Note: See TracChangeset for help on using the changeset viewer.