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 trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90 @ 4765

Last change on this file since 4765 was 4333, checked in by clem, 10 years ago

remove remaining bugs in LIM3, so that it can run in both regional and global config
  • Property svn:keywords set to Id
File size: 8.2 KB
Line 
1MODULE limhdf
2   !!======================================================================
3   !!                    ***  MODULE limhdf   ***
4   !! LIM ice model : horizontal diffusion of sea-ice quantities
5   !!======================================================================
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   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                      LIM3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_hdf       : diffusion trend on sea-ice variable
15   !!----------------------------------------------------------------------
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) 
24
25   IMPLICIT NONE
26   PRIVATE
27
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
32   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   efact   ! metric coefficient
33
34   !! * Substitution
35#  include "vectopt_loop_substitute.h90"
36   !!----------------------------------------------------------------------
37   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
38   !! $Id$
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE lim_hdf( ptab )
44      !!-------------------------------------------------------------------
45      !!                  ***  ROUTINE lim_hdf  ***
46      !!
47      !! ** purpose :   Compute and add the diffusive trend on sea-ice variables
48      !!
49      !! ** method  :   Second order diffusive operator evaluated using a
50      !!              Cranck-Nicholson time Scheme.
51      !!
52      !! ** Action  :    update ptab with the diffusive contribution
53      !!-------------------------------------------------------------------
54      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   ptab    ! Field on which the diffusion is applied
55      !
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
61      !!-------------------------------------------------------------------
62     
63      CALL wrk_alloc( jpi, jpj, zrlx, zflu, zflv, zdiv0, zdiv, ztab0 )
64
65      !                       !==  Initialisation  ==!
66      !
67      IF( linit ) THEN              ! Metric coefficient (compute at the first call and saved in efact)
68         ALLOCATE( efact(jpi,jpj) , STAT=ierr )
69         IF( lk_mpp    )   CALL mpp_sum( ierr )
70         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'lim_hdf : unable to allocate arrays' )
71         DO jj = 2, jpjm1 
72            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) )
74            END DO
75         END DO
76         linit = .FALSE.
77      ENDIF
78      !                             ! 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
81      !
82      ztab0(:, : ) = ptab(:,:)      ! Arrays initialization
83      zdiv0(:, 1 ) = 0._wp
84      zdiv0(:,jpj) = 0._wp
85      IF( .NOT.lk_vopt_loop ) THEN
86         zflu (jpi,:) = 0._wp   
87         zflv (jpi,:) = 0._wp
88         zdiv0(1,  :) = 0._wp
89         zdiv0(jpi,:) = 0._wp
90      ENDIF
91
92      zconv = 1._wp           !==  horizontal diffusion using a Crant-Nicholson scheme  ==!
93      iter  = 0
94      !
95      DO WHILE( zconv > ( 2._wp * epsi04 ) .AND. iter <= its )   ! Sub-time step loop
96         !
97         iter = iter + 1                                 ! incrementation of the sub-time step number
98         !
99         DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
100            DO ji = 1 , fs_jpim1   ! vector opt.
101               zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
102               zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) / e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )
103            END DO
104         END DO
105         !
106         DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes
107            DO ji = fs_2 , fs_jpim1   ! vector opt.
108               zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   &
109                  &            + zflv(ji,jj) - zflv(ji  ,jj-1)  ) / ( e1t (ji,jj) * e2t (ji,jj) )
110            END DO
111         END DO
112         !
113         IF( iter == 1 )   zdiv0(:,:) = zdiv(:,:)        ! save the 1st evaluation of the diffusive trend in zdiv0
114         !
115         DO jj = 2, jpjm1                                ! iterative evaluation
116            DO ji = fs_2 , fs_jpim1   ! vector opt.
117               zrlxint = (   ztab0(ji,jj)    &
118                  &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + efact(ji,jj) * 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) + om * ( zrlxint - ptab(ji,jj) )
122            END DO
123         END DO
124         CALL lbc_lnk( zrlx, 'T', 1. )                   ! lateral boundary condition
125         !
126         zconv = 0._wp                                   ! convergence test
127         DO jj = 2, jpjm1
128            DO ji = fs_2, fs_jpim1
129               zconv = MAX( zconv, ABS( zrlx(ji,jj) - ptab(ji,jj) )  )
130            END DO
131         END DO
132         IF( lk_mpp )   CALL mpp_max( zconv )            ! max over the global domain
133         !
134         ptab(:,:) = zrlx(:,:)
135         !
136      END DO                                       ! end of sub-time step loop
137
138      ! -----------------------
139      !!! final step (clem) !!!
140      DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
141         DO ji = 1 , fs_jpim1   ! vector opt.
142            zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
143            zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) / e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )
144         END DO
145      END DO
146      !
147      DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes
148         DO ji = fs_2 , fs_jpim1   ! vector opt.
149            zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   &
150                 &            + zflv(ji,jj) - zflv(ji  ,jj-1)  ) / ( e1t (ji,jj) * e2t (ji,jj) )
151            ptab(ji,jj) = ztab0(ji,jj) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj) )
152         END DO
153      END DO
154      CALL lbc_lnk( ptab, 'T', 1. )                   ! lateral boundary condition
155      !!! final step (clem) !!!
156      ! -----------------------
157
158      IF(ln_ctl)   THEN
159         zrlx(:,:) = ptab(:,:) - ztab0(:,:)
160         WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter
161         CALL prt_ctl( tab2d_1=zrlx, clinfo1=charout )
162      ENDIF
163      !
164      CALL wrk_dealloc( jpi, jpj, zrlx, zflu, zflv, zdiv0, zdiv, ztab0 )
165      !
166   END SUBROUTINE lim_hdf
167
168#else
169   !!----------------------------------------------------------------------
170   !!   Default option          Dummy module           NO LIM sea-ice model
171   !!----------------------------------------------------------------------
172CONTAINS
173   SUBROUTINE lim_hdf         ! Empty routine
174   END SUBROUTINE lim_hdf
175#endif
176
177   !!======================================================================
178END MODULE limhdf
Note: See TracBrowser for help on using the repository browser.