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/NEMO/LIM_SRC – NEMO

source: trunk/NEMO/LIM_SRC/limhdf.F90 @ 12

Last change on this file since 12 was 12, checked in by opalod, 21 years ago

CT : UPDATE001 : First major NEMO update
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.6 KB
Line 
1MODULE limhdf
2#if defined key_ice_lim
3   !!======================================================================
4   !!                    ***  MODULE limhdf   ***
5   !! LIM diffusion ice model : sea-ice variables horizontal diffusion
6   !!======================================================================
7
8   !!----------------------------------------------------------------------
9   !!   lim_hdf  : diffusion trend on sea-ice variable
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE dom_oce
13   USE ice_oce         ! ice variables
14   USE in_out_manager
15   USE ice
16!  USE limdyn
17   USE lbclnk
18   USE lib_mpp
19
20   IMPLICIT NONE
21   PRIVATE
22
23   !! * Routine accessibility
24   PUBLIC lim_hdf    ! called by lim_tra
25
26   !! * Module variables
27   LOGICAL  ::   linit = .TRUE.              ! ???
28   REAL(wp) ::   epsi04 = 1e-04              ! constant
29   REAL(wp), DIMENSION(jpi,jpj) ::   zfact   ! ???
30
31   !! * Substitution
32#  include "vectopt_loop_substitute.h90"
33   !!----------------------------------------------------------------------
34   !!   LIM 2.0 , UCL-LODYC-IPSL  (2003)
35   !!----------------------------------------------------------------------
36
37CONTAINS
38
39   SUBROUTINE lim_hdf( ptab )
40      !!-------------------------------------------------------------------
41      !!                  ***  ROUTINE lim_hdf  ***
42      !!
43      !! ** purpose :   Compute and add the diffusive trend on sea-ice
44      !!      variables
45      !!
46      !! ** method  :   Second order diffusive operator evaluated using a
47      !!      Cranck-Nicholson time Scheme.
48      !!
49      !! ** Action  :    update ptab with the diffusive contribution
50      !!
51      !! History :
52      !!        !  00-01 (LIM) Original code
53      !!        !  01-05 (G. Madec, R. Hordoir) opa norm
54      !!        !  02-08 (C. Ethe)  F90, free form
55      !!-------------------------------------------------------------------
56      ! * Arguments
57      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   &
58         ptab                 ! Field on which the diffusion is applied 
59      REAL(wp), DIMENSION(jpi,jpj) ::   &
60         ptab0                ! ???
61
62      ! * Local variables
63      INTEGER ::  ji, jj      ! dummy loop indices
64      INTEGER ::  &
65         its, iter            ! temporary integers
66      REAL(wp) ::  &
67         zalfa, zrlxint, zconv, zeps   ! temporary scalars
68      REAL(wp), DIMENSION(jpi,jpj) ::  & 
69         zrlx, zflu, zflv, &  ! temporary workspaces
70         zdiv0, zdiv          !    "           "
71      !!-------------------------------------------------------------------
72
73      ! Initialisation
74      ! ---------------   
75      ! Time integration parameters
76      zalfa = 0.5       ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit
77      its   = 100       ! Maximum number of iteration
78      zeps  =  2. * epsi04
79
80      ! Arrays initialization
81      ptab0 (:, : ) = ptab(:,:)
82!bug  zflu (:,jpj) = 0.e0
83!bug  zflv (:,jpj) = 0.e0
84      zdiv0(:, 1 ) = 0.e0
85      zdiv0(:,jpj) = 0.e0
86      IF( .NOT.lk_vopt_loop ) THEN
87         zflu (jpi,:) = 0.e0   
88         zflv (jpi,:) = 0.e0
89         zdiv0(1,  :) = 0.e0
90         zdiv0(jpi,:) = 0.e0
91      ENDIF
92
93      ! Metric coefficient (compute at the first call and saved in
94      IF( linit ) THEN
95         DO jj = 2, jpjm1 
96            DO ji = fs_2 , fs_jpim1   ! vector opt.
97               zfact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj  ) + e1v(ji,jj) + e1v(ji,jj-1) ) &
98                             / ( e1t(ji,jj) * e2t(ji,jj) )
99            END DO
100         END DO
101         linit = .FALSE.
102      ENDIF
103
104
105      ! Sub-time step loop
106      zconv = 1.e0
107      iter  = 0
108
109      !                                                   !===================
110      DO WHILE ( ( zconv > zeps ) .AND. (iter <= its) )   ! Sub-time step loop
111         !                                                !===================
112         ! incrementation of the sub-time step number
113         iter = iter + 1
114
115         ! diffusive fluxes in U- and V- direction
116         DO jj = 1, jpjm1
117            DO ji = 1 , fs_jpim1   ! vector opt.
118               zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
119               zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) / e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )
120            END DO
121         END DO
122
123         ! diffusive trend : divergence of the fluxes
124         DO jj= 2, jpjm1
125            DO ji = fs_2 , fs_jpim1   ! vector opt.
126               zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   &
127                  &            + zflv(ji,jj) - zflv(ji  ,jj-1)  ) / ( e1t (ji,jj) * e2t (ji,jj) )
128            END DO
129         END DO
130
131         ! save the first evaluation of the diffusive trend in zdiv0
132         IF( iter == 1 )   zdiv0(:,:) = zdiv(:,:)       
133
134         ! XXXX iterative evaluation?????
135         DO jj = 2, jpjm1
136            DO ji = fs_2 , fs_jpim1   ! vector opt.
137               zrlxint = (   ptab0(ji,jj)    &
138                  &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + zfact(ji,jj) * ptab(ji,jj) )   &
139                  &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )  )                             & 
140                  &    / ( 1.0 + zalfa * rdt_ice * zfact(ji,jj) )
141               zrlx(ji,jj) = ptab(ji,jj) + om * ( zrlxint - ptab(ji,jj) )
142            END DO
143         END DO
144
145         ! convergence test
146         zconv = 0.0
147         DO jj = 2, jpjm1
148            DO ji = 2, jpim1
149               zconv = MAX( zconv, ABS( zrlx(ji,jj) - ptab(ji,jj) )  )
150            END DO
151         END DO
152         IF( lk_mpp )   CALL mpp_max( zconv )   ! max over the global domain
153
154         DO jj = 2, jpjm1
155            DO ji = 2 , jpim1
156               ptab(ji,jj) = zrlx(ji,jj)
157            END DO
158         END DO
159
160         ! lateral boundary condition on ptab
161         CALL lbc_lnk( ptab, 'T', 1. )
162         !                                         !==========================
163      END DO                                       ! end of sub-time step loop
164      !                                            !==========================
165
166      ptab(:,:) = ptab(:,:)
167      IF( l_ctl .AND. lwp ) THEN
168         WRITE(numout,*) ' lim_hdf  : ', SUM( ptab-ptab0 ), ' zconv= ', zconv, ' iter= ', iter
169      ENDIF
170
171   END SUBROUTINE lim_hdf
172#else
173   !!======================================================================
174   !!                       ***  MODULE limhdf   ***
175   !!                          no sea ice model
176   !!======================================================================
177CONTAINS
178   SUBROUTINE lim_hdf         ! Empty routine
179   END SUBROUTINE lim_hdf
180#endif
181
182   !!======================================================================
183END MODULE limhdf
Note: See TracBrowser for help on using the repository browser.