source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90 @ 4333

Last change on this file since 4333 was 4333, checked in by clem, 7 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.