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 @ 2715

Last change on this file since 2715 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk

  • Property svn:keywords set to Id
File size: 7.4 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 prtctl           ! Print control
21   USE in_out_manager   ! I/O manager
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   lim_hdf     ! called by lim_tra
27
28   LOGICAL  ::   linit = .TRUE.              ! initialization flag (set to flase after the 1st call)
29   REAL(wp) ::   epsi04 = 1e-04              ! constant
30   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   efact   ! metric coefficient
31
32   !! * Substitution
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
36   !! $Id$
37   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE lim_hdf( ptab )
42      !!-------------------------------------------------------------------
43      !!                  ***  ROUTINE lim_hdf  ***
44      !!
45      !! ** purpose :   Compute and add the diffusive trend on sea-ice variables
46      !!
47      !! ** method  :   Second order diffusive operator evaluated using a
48      !!              Cranck-Nicholson time Scheme.
49      !!
50      !! ** Action  :    update ptab with the diffusive contribution
51      !!-------------------------------------------------------------------
52      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
53      USE wrk_nemo, ONLY:   zflu => wrk_2d_11, zdiv  => wrk_2d_13, zrlx  => wrk_2d_15 
54      USE wrk_nemo, ONLY:   zflv => wrk_2d_12, zdiv0 => wrk_2d_14, ztab0 => wrk_2d_16
55      !
56      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   ptab    ! Field on which the diffusion is applied
57      !
58      INTEGER  ::  ji, jj            ! dummy loop indices
59      INTEGER  ::  its, iter, ierr   ! local integers
60      REAL(wp) ::   zalfa, zrlxint, zconv, zeps   ! local scalars
61      CHARACTER(lc) ::   charout   ! local character
62      !!-------------------------------------------------------------------
63     
64      IF( wrk_in_use(2, 11,12,13,14,15,16) ) THEN
65         CALL ctl_stop( 'lim_hdf: requested workspace arrays unavailable' )   ;   RETURN
66      ENDIF
67
68      !                       !==  Initialisation  ==!
69      !
70      IF( linit ) THEN              ! Metric coefficient (compute at the first call and saved in efact)
71         ALLOCATE( efact(jpi,jpj) , STAT=ierr )
72         IF( lk_mpp    )   CALL mpp_sum( ierr )
73         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'lim_hdf : unable to allocate arrays' )
74         DO jj = 2, jpjm1 
75            DO ji = fs_2 , fs_jpim1   ! vector opt.
76               efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) )
77            END DO
78         END DO
79         linit = .FALSE.
80      ENDIF
81      !                             ! Time integration parameters
82      zalfa = 0.5_wp                      ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit
83      its   = 100                         ! Maximum number of iteration
84      zeps  =  2._wp * epsi04
85      !
86      ztab0(:, : ) = ptab(:,:)      ! Arrays initialization
87      zdiv0(:, 1 ) = 0._wp
88      zdiv0(:,jpj) = 0._wp
89      IF( .NOT.lk_vopt_loop ) THEN
90         zflu (jpi,:) = 0._wp   
91         zflv (jpi,:) = 0._wp
92         zdiv0(1,  :) = 0._wp
93         zdiv0(jpi,:) = 0._wp
94      ENDIF
95
96      zconv = 1._wp           !==  horizontal diffusion using a Crant-Nicholson scheme  ==!
97      iter  = 0
98      !
99      DO WHILE( zconv > zeps .AND. iter <= its )   ! Sub-time step loop
100         !
101         iter = iter + 1                                 ! incrementation of the sub-time step number
102         !
103         DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
104            DO ji = 1 , fs_jpim1   ! vector opt.
105               zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
106               zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) / e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )
107            END DO
108         END DO
109         !
110         DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes
111            DO ji = fs_2 , fs_jpim1   ! vector opt.
112               zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   &
113                  &            + zflv(ji,jj) - zflv(ji  ,jj-1)  ) / ( e1t (ji,jj) * e2t (ji,jj) )
114            END DO
115         END DO
116         !
117         IF( iter == 1 )   zdiv0(:,:) = zdiv(:,:)        ! save the 1st evaluation of the diffusive trend in zdiv0
118         !
119         DO jj = 2, jpjm1                                ! iterative evaluation
120            DO ji = fs_2 , fs_jpim1   ! vector opt.
121               zrlxint = (   ztab0(ji,jj)    &
122                  &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj) )   &
123                  &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )  )                             & 
124                  &    / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) )
125               zrlx(ji,jj) = ptab(ji,jj) + om * ( zrlxint - ptab(ji,jj) )
126            END DO
127         END DO
128         CALL lbc_lnk( zrlx, 'T', 1. )                   ! lateral boundary condition
129         !
130         zconv = 0._wp                                   ! convergence test
131         DO jj = 2, jpjm1
132            DO ji = fs_2, fs_jpim1
133               zconv = MAX( zconv, ABS( zrlx(ji,jj) - ptab(ji,jj) )  )
134            END DO
135         END DO
136         IF( lk_mpp )   CALL mpp_max( zconv )            ! max over the global domain
137         !
138         ptab(:,:) = zrlx(:,:)
139         !
140      END DO                                       ! end of sub-time step loop
141
142      IF(ln_ctl)   THEN
143         zrlx(:,:) = ptab(:,:) - ztab0(:,:)
144         WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter
145         CALL prt_ctl( tab2d_1=zrlx, clinfo1=charout )
146      ENDIF
147      !
148      IF( wrk_not_released(2, 11,12,13,14,15,16) )   CALL ctl_stop('lim_hdf: failed to release workspace arrays')
149      !
150   END SUBROUTINE lim_hdf
151
152#else
153   !!----------------------------------------------------------------------
154   !!   Default option          Dummy module           NO LIM sea-ice model
155   !!----------------------------------------------------------------------
156CONTAINS
157   SUBROUTINE lim_hdf         ! Empty routine
158   END SUBROUTINE lim_hdf
159#endif
160
161   !!======================================================================
162END MODULE limhdf
Note: See TracBrowser for help on using the repository browser.