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.
limhdf.F90 in branches/2015/dev_r5302_CNRS18_HPC_scalability/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2015/dev_r5302_CNRS18_HPC_scalability/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90 @ 5320

Last change on this file since 5320 was 5320, checked in by mcastril, 9 years ago

ticket #1523 Convergence Check
  • Property svn:keywords set to Id
File size: 9.1 KB
RevLine 
[825]1MODULE limhdf
2   !!======================================================================
3   !!                    ***  MODULE limhdf   ***
4   !! LIM ice model : horizontal diffusion of sea-ice quantities
5   !!======================================================================
[2715]6   !! History :  LIM  !  2000-01 (LIM) Original code
7   !!             -   !  2001-05 (G. Madec, R. Hordoir) opa norm
8   !!            1.0  !  2002-08 (C. Ethe)  F90, free form
9   !!----------------------------------------------------------------------
[825]10#if defined key_lim3
11   !!----------------------------------------------------------------------
[834]12   !!   'key_lim3'                                      LIM3 sea-ice model
[825]13   !!----------------------------------------------------------------------
[3625]14   !!   lim_hdf       : diffusion trend on sea-ice variable
[825]15   !!----------------------------------------------------------------------
[3625]16   USE dom_oce        ! ocean domain
17   USE ice            ! LIM-3: ice variables
18   USE lbclnk         ! lateral boundary condition - MPP exchanges
19   USE lib_mpp        ! MPP library
20   USE wrk_nemo       ! work arrays
21   USE prtctl         ! Print control
22   USE in_out_manager ! I/O manager
23   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
[825]24
25   IMPLICIT NONE
26   PRIVATE
27
[5123]28   PUBLIC   lim_hdf     ! called by lim_trp
[825]29
[5123]30   LOGICAL  ::   linit = .TRUE.                             ! initialization flag (set to flase after the 1st call)
[2715]31   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   efact   ! metric coefficient
[825]32
33   !! * Substitution
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
[4161]36   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
[1156]37   !! $Id$
[2715]38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[825]39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE lim_hdf( ptab )
43      !!-------------------------------------------------------------------
44      !!                  ***  ROUTINE lim_hdf  ***
45      !!
[2715]46      !! ** purpose :   Compute and add the diffusive trend on sea-ice variables
[825]47      !!
48      !! ** method  :   Second order diffusive operator evaluated using a
[2715]49      !!              Cranck-Nicholson time Scheme.
[825]50      !!
51      !! ** Action  :    update ptab with the diffusive contribution
52      !!-------------------------------------------------------------------
[2715]53      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   ptab    ! Field on which the diffusion is applied
54      !
[5123]55      INTEGER                           ::  ji, jj                    ! dummy loop indices
56      INTEGER                           ::  iter, ierr           ! local integers
57      REAL(wp)                          ::  zrlxint, zconv     ! local scalars
58      REAL(wp), POINTER, DIMENSION(:,:) ::  zrlx, zflu, zflv, zdiv0, zdiv, ztab0
59      CHARACTER(lc)                     ::  charout                   ! local character
60      REAL(wp), PARAMETER               ::  zrelax = 0.5_wp           ! relaxation constant for iterative procedure
61      REAL(wp), PARAMETER               ::  zalfa  = 0.5_wp           ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit
62      INTEGER , PARAMETER               ::  its    = 100              ! Maximum number of iteration
[825]63      !!-------------------------------------------------------------------
[2715]64     
[3294]65      CALL wrk_alloc( jpi, jpj, zrlx, zflu, zflv, zdiv0, zdiv, ztab0 )
[825]66
[2715]67      !                       !==  Initialisation  ==!
68      !
69      IF( linit ) THEN              ! Metric coefficient (compute at the first call and saved in efact)
70         ALLOCATE( efact(jpi,jpj) , STAT=ierr )
71         IF( lk_mpp    )   CALL mpp_sum( ierr )
72         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'lim_hdf : unable to allocate arrays' )
[825]73         DO jj = 2, jpjm1 
74            DO ji = fs_2 , fs_jpim1   ! vector opt.
[5123]75               efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) * r1_e12t(ji,jj)
[825]76            END DO
77         END DO
78         linit = .FALSE.
79      ENDIF
[2715]80      !                             ! Time integration parameters
81      !
82      ztab0(:, : ) = ptab(:,:)      ! Arrays initialization
83      zdiv0(:, 1 ) = 0._wp
84      zdiv0(:,jpj) = 0._wp
[4990]85      zflu (jpi,:) = 0._wp   
86      zflv (jpi,:) = 0._wp
87      zdiv0(1,  :) = 0._wp
88      zdiv0(jpi,:) = 0._wp
[825]89
[2715]90      zconv = 1._wp           !==  horizontal diffusion using a Crant-Nicholson scheme  ==!
[825]91      iter  = 0
[2715]92      !
[5123]93      DO WHILE( zconv > ( 2._wp * 1.e-04 ) .AND. iter <= its )   ! Sub-time step loop
[2715]94         !
95         iter = iter + 1                                 ! incrementation of the sub-time step number
96         !
97         DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
[825]98            DO ji = 1 , fs_jpim1   ! vector opt.
[5123]99               zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
100               zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )
[825]101            END DO
102         END DO
[2715]103         !
104         DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes
[825]105            DO ji = fs_2 , fs_jpim1   ! vector opt.
[5123]106               zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj)
[825]107            END DO
108         END DO
[2715]109         !
110         IF( iter == 1 )   zdiv0(:,:) = zdiv(:,:)        ! save the 1st evaluation of the diffusive trend in zdiv0
111         !
112         DO jj = 2, jpjm1                                ! iterative evaluation
[825]113            DO ji = fs_2 , fs_jpim1   ! vector opt.
[2715]114               zrlxint = (   ztab0(ji,jj)    &
115                  &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj) )   &
[5123]116                  &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )                               & 
117                  &      ) / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) )
118               zrlx(ji,jj) = ptab(ji,jj) + zrelax * ( zrlxint - ptab(ji,jj) )
[825]119            END DO
120         END DO
[2715]121         CALL lbc_lnk( zrlx, 'T', 1. )                   ! lateral boundary condition
122         !
[5320]123
124         !! The convergence check optimization can be activated in the namelist using
125         !! the parameters cc_opt to activate the optimization and cc_freq to set the frequency of
126         !! Convergence Checks. These parameters must go into namicedyn section.
127
128         IF ( cc_opt == .true. ) THEN
129            !Convergence test every cc_freq time-steps
130            IF ( MOD( iter-1 , cc_freq ) == 0 )  THEN
131               zconv = 0._wp                                   ! convergence test
132               DO jj = 2, jpjm1
133                  DO ji = fs_2, fs_jpim1
134                     zconv = MAX( zconv, ABS( zrlx(ji,jj) - ptab(ji,jj) )  )
135                  END DO
136               END DO
137               IF( lk_mpp )   CALL mpp_max( zconv )            ! max over the global domain             
138            ENDIF
139         ELSE
140             zconv = 0._wp                                   ! convergence test
141             DO jj = 2, jpjm1
142                DO ji = fs_2, fs_jpim1
143                   zconv = MAX( zconv, ABS( zrlx(ji,jj) - ptab(ji,jj) )  )
144                END DO
145             END DO
146             IF( lk_mpp )   CALL mpp_max( zconv )            ! max over the global domain
147         END IF
148         
[2715]149         !
150         ptab(:,:) = zrlx(:,:)
151         !
[825]152      END DO                                       ! end of sub-time step loop
153
[4161]154      ! -----------------------
155      !!! final step (clem) !!!
156      DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
157         DO ji = 1 , fs_jpim1   ! vector opt.
[5123]158            zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
159            zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )
[4161]160         END DO
161      END DO
162      !
163      DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes
164         DO ji = fs_2 , fs_jpim1   ! vector opt.
[5123]165            zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj)
[4161]166            ptab(ji,jj) = ztab0(ji,jj) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj) )
167         END DO
168      END DO
169      CALL lbc_lnk( ptab, 'T', 1. )                   ! lateral boundary condition
170      !!! final step (clem) !!!
171      ! -----------------------
172
[825]173      IF(ln_ctl)   THEN
[2715]174         zrlx(:,:) = ptab(:,:) - ztab0(:,:)
[825]175         WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter
[2715]176         CALL prt_ctl( tab2d_1=zrlx, clinfo1=charout )
[825]177      ENDIF
[2715]178      !
[3294]179      CALL wrk_dealloc( jpi, jpj, zrlx, zflu, zflv, zdiv0, zdiv, ztab0 )
[2715]180      !
[825]181   END SUBROUTINE lim_hdf
182
183#else
184   !!----------------------------------------------------------------------
185   !!   Default option          Dummy module           NO LIM sea-ice model
186   !!----------------------------------------------------------------------
187CONTAINS
188   SUBROUTINE lim_hdf         ! Empty routine
189   END SUBROUTINE lim_hdf
190#endif
191
192   !!======================================================================
193END MODULE limhdf
Note: See TracBrowser for help on using the repository browser.