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/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90 @ 7525

Last change on this file since 7525 was 7525, checked in by mocavero, 7 years ago

changes after review
  • Property svn:keywords set to Id
File size: 8.2 KB
Line 
1MODULE limhdf_2
2   !!======================================================================
3   !!                    ***  MODULE limhdf_2   ***
4   !! LIM 2.0 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_lim2
11   !!----------------------------------------------------------------------
12   !!   'key_lim2'                                    LIM 2.0 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_hdf_2  : diffusion trend on sea-ice variable
15   !!----------------------------------------------------------------------
16   USE dom_oce          ! ocean domain
17   USE ice_2            ! LIM-2: 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_2         ! called by limtrp_2.F90
29
30   LOGICAL  ::   linit = .TRUE.   ! ! initialization flag (set to flase after the 1st call)
31   REAL(wp) ::   epsi04 = 1e-04   ! constant
32   
33   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   efact   ! metric coefficient
34
35   !! * Substitution
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/LIM2 4.0 , UCL - NEMO Consortium (2010)
39   !! $Id$
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE lim_hdf_2( ptab )
45      !!-------------------------------------------------------------------
46      !!                  ***  ROUTINE lim_hdf_2  ***
47      !!
48      !! ** purpose :   Compute and add the diffusive trend on sea-ice variables
49      !!
50      !! ** method  :   Second order diffusive operator evaluated using a
51      !!              Cranck-Nicholson time Scheme.
52      !!
53      !! ** Action  :    update ptab with the diffusive contribution
54      !!-------------------------------------------------------------------
55      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   ptab   ! Field on which the diffusion is applied 
56      !
57      INTEGER  ::   ji, jj            ! dummy loop indices
58      INTEGER  ::   its, iter, ierr   ! local integers
59      REAL(wp) ::   zalfa, zrlxint, zconv, zeps   ! local scalars
60      REAL(wp), DIMENSION(:,:), POINTER :: zrlx, zflu, zflv, zdiv0, zdiv, ztab0 
61      CHARACTER (len=55) :: charout
62      !!-------------------------------------------------------------------
63
64      CALL wrk_alloc( jpi, jpj, zrlx, zflu, zflv, zdiv0, zdiv, ztab0 )
65
66      !                       !==  Initialisation  ==!
67      !
68      IF( linit ) THEN              ! Metric coefficient (compute at the first call and saved in efact)
69         ALLOCATE( efact(jpi,jpj) , STAT=ierr )
70         IF( lk_mpp    )   CALL mpp_sum( ierr )
71         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'lim_hdf_2 : unable to allocate standard arrays' )
72!$OMP PARALLEL DO schedule(static) private(jj, ji)
73         DO jj = 2, jpjm1 
74            DO ji = fs_2 , fs_jpim1   ! vector opt.
75               efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) )
76            END DO
77         END DO
78         linit = .FALSE.
79      ENDIF
80      !
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!$OMP PARALLEL DO schedule(static) private(jj, ji)
87      DO jj = 1, jpj
88         DO ji = 1, jpi
89            ztab0(ji, jj ) = ptab(ji,jj)      ! Arrays initialization
90         END DO
91      END DO
92      zdiv0(:, 1 ) = 0._wp
93      zdiv0(:,jpj) = 0._wp
94      zflu (jpi,:) = 0._wp
95      zflv (jpi,:) = 0._wp
96      zdiv0(1,  :) = 0._wp
97      zdiv0(jpi,:) = 0._wp
98
99      zconv = 1._wp           !==  horizontal diffusion using a Crant-Nicholson scheme  ==!
100      iter  = 0
101      !
102      DO WHILE (  zconv > zeps  .AND.  iter <= its  )    ! Sub-time step loop
103         !
104         iter = iter + 1                                       ! incrementation of the sub-time step number
105         !
106!$OMP PARALLEL
107!$OMP DO schedule(static) private(jj, ji)
108         DO jj = 1, jpjm1                                      ! diffusive fluxes in U- and V- direction
109            DO ji = 1 , fs_jpim1   ! vector opt.
110               zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
111               zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) / e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )
112            END DO
113         END DO
114         !
115!$OMP DO schedule(static) private(jj, ji)
116         DO jj= 2, jpjm1                                       ! diffusive trend : divergence of the fluxes
117            DO ji = fs_2 , fs_jpim1   ! vector opt.
118               zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   &
119                  &            + zflv(ji,jj) - zflv(ji  ,jj-1)  ) / ( e1t (ji,jj) * e2t (ji,jj) )
120            END DO
121         END DO
122!$OMP END DO NOWAIT
123!$OMP END PARALLEL
124         !
125         IF( iter == 1 ) THEN
126!$OMP PARALLEL DO schedule(static) private(jj, ji)
127             DO jj= 2, jpjm1
128                DO ji = fs_2 , fs_jpim1   ! vector opt.
129                   zdiv0(ji,jj) = zdiv(ji,jj)              ! save the 1st evaluation of the diffusive trend in zdiv0
130                END DO
131             END DO
132         END IF
133         !
134!$OMP PARALLEL DO schedule(static) private(jj,ji,zrlxint)
135         DO jj = 2, jpjm1                                      ! iterative evaluation
136            DO ji = fs_2 , fs_jpim1   ! vector opt.
137               zrlxint = (   ztab0(ji,jj)    &
138                  &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj) )   &
139                  &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )  )                             & 
140                  &    / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) )
141               zrlx(ji,jj) = ptab(ji,jj) + om * ( zrlxint - ptab(ji,jj) )
142            END DO
143         END DO
144         CALL lbc_lnk( zrlx, 'T', 1. )                         ! lateral boundary condition
145
146         zconv = 0._wp                                         ! convergence test
147
148!$OMP PARALLEL DO schedule(static) private(jj,ji) REDUCTION(MAX:zconv)
149         DO jj = 2, jpjm1
150            DO ji = 2, jpim1
151               zconv = MAX( zconv, ABS( zrlx(ji,jj) - ptab(ji,jj) )  )
152            END DO
153         END DO
154         IF( lk_mpp )   CALL mpp_max( zconv )                  ! max over the global domain
155
156!$OMP PARALLEL DO schedule(static) private(jj, ji)
157         DO jj = 1, jpj
158            DO ji = 1, jpi
159               ptab(ji,jj) = zrlx(ji,jj)
160            END DO
161         END DO
162         !
163      END DO                                             ! end of sub-time step loop
164
165      IF(ln_ctl)   THEN
166!$OMP PARALLEL DO schedule(static) private(jj, ji)
167         DO jj = 1, jpj
168            DO ji = 1, jpi
169               zrlx(ji,jj) = ptab(ji,jj) - ztab0(ji,jj)
170            END DO
171         END DO
172         WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter
173         CALL prt_ctl( tab2d_1=zrlx, clinfo1=charout )
174      ENDIF
175      !
176      CALL wrk_dealloc( jpi, jpj, zrlx, zflu, zflv, zdiv0, zdiv, ztab0 )
177      !
178   END SUBROUTINE lim_hdf_2
179
180#else
181   !!----------------------------------------------------------------------
182   !!   Default option          Dummy module       NO LIM 2.0 sea-ice model
183   !!----------------------------------------------------------------------
184CONTAINS
185   SUBROUTINE lim_hdf_2       ! Empty routine
186   END SUBROUTINE lim_hdf_2
187#endif
188
189   !!======================================================================
190END MODULE limhdf_2
Note: See TracBrowser for help on using the repository browser.