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_2.F90 in branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 13 years ago

update licence of all NEMO files...

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