source: branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90

Last change on this file was 7993, checked in by frrh, 4 years ago

Merge in missing revisions 6428:2477 inclusive and 6482 from nemo_v3_6_STABLE
branch. In ptic, this includes the fix for restartability of runoff fields in coupled
models. Evolution of coupled models will therefor be affected.

These changes donot affect evolution of the current stand-alone NEMO-CICE GO6
standard configuration.

Work and testing documented in Met Office GMED ticket 320.

File size: 12.6 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   !!            3.0  !  2015-08 (O. Tintó and M. Castrillo)  added lim_hdf (multiple)
10   !!----------------------------------------------------------------------
11#if defined key_lim3
12   !!----------------------------------------------------------------------
13   !!   'key_lim3'                                      LIM3 sea-ice model
14   !!----------------------------------------------------------------------
15   !!   lim_hdf       : diffusion trend on sea-ice variable
16   !!   lim_hdf_init  : initialisation of diffusion trend on sea-ice variable
17   !!----------------------------------------------------------------------
18   USE dom_oce        ! ocean domain
19   USE ice            ! LIM-3: ice variables
20   USE lbclnk         ! lateral boundary condition - MPP exchanges
21   USE lib_mpp        ! MPP library
22   USE wrk_nemo       ! work arrays
23   USE prtctl         ! Print control
24   USE in_out_manager ! I/O manager
25   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   lim_hdf ! called by lim_trp
31   PUBLIC   lim_hdf_init    ! called by sbc_lim_init
32
33   LOGICAL  ::   linit = .TRUE.                             ! initialization flag (set to flase after the 1st call)
34   INTEGER  ::   nn_convfrq                                 !:  convergence check frequency of the Crant-Nicholson scheme
35   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   efact   ! metric coefficient
36
37   !! * Substitution
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE lim_hdf( ptab , ihdf_vars , jpl , nlay_i )
47      !!-------------------------------------------------------------------
48      !!                  ***  ROUTINE lim_hdf  ***
49      !!
50      !! ** purpose :   Compute and add the diffusive trend on sea-ice variables
51      !!
52      !! ** method  :   Second order diffusive operator evaluated using a
53      !!              Cranck-Nicholson time Scheme.
54      !!
55      !! ** Action  :    update ptab with the diffusive contribution
56      !!-------------------------------------------------------------------
57      INTEGER                           :: jpl, nlay_i, isize, ihdf_vars
58      REAL(wp),  DIMENSION(:,:,:), INTENT( inout ),TARGET ::   ptab    ! Field on which the diffusion is applied
59      !
60      INTEGER                           ::  ji, jj, jk, jl , jm               ! dummy loop indices
61      INTEGER                           ::  iter, ierr           ! local integers
62      REAL(wp)                          ::  zrlxint     ! local scalars
63      REAL(wp), POINTER , DIMENSION ( : )        :: zconv     ! local scalars
64      REAL(wp), POINTER , DIMENSION(:,:,:) ::  zrlx,zdiv0, ztab0
65      REAL(wp), POINTER , DIMENSION(:,:) ::  zflu, zflv, zdiv
66      CHARACTER(lc)                     ::  charout                   ! local character
67      REAL(wp), PARAMETER               ::  zrelax = 0.5_wp           ! relaxation constant for iterative procedure
68      REAL(wp), PARAMETER               ::  zalfa  = 0.5_wp           ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit
69      INTEGER , PARAMETER               ::  its    = 100              ! Maximum number of iteration
70      !!-------------------------------------------------------------------
71      TYPE(arrayptr)   , ALLOCATABLE, DIMENSION(:) ::   pt2d_array, zrlx_array
72      CHARACTER(len=1) , ALLOCATABLE, DIMENSION(:) ::   type_array ! define the nature of ptab array grid-points
73      !                                                            ! = T , U , V , F , W and I points
74      REAL(wp)        , ALLOCATABLE, DIMENSION(:)  ::   psgn_array    ! =-1 the sign change across the north fold boundary
75
76     !!---------------------------------------------------------------------
77
78      !                       !==  Initialisation  ==!
79      ! +1 open water diffusion
80      isize = jpl*(ihdf_vars+nlay_i)+1
81      ALLOCATE( zconv (isize) )
82      ALLOCATE( pt2d_array(isize) , zrlx_array(isize) )
83      ALLOCATE( type_array(isize) )
84      ALLOCATE( psgn_array(isize) )
85     
86      CALL wrk_alloc( jpi, jpj, isize, zrlx, zdiv0, ztab0 )
87      CALL wrk_alloc( jpi, jpj, zflu, zflv, zdiv )
88
89      DO jk= 1 , isize
90         pt2d_array(jk)%pt2d=>ptab(:,:,jk)
91         zrlx_array(jk)%pt2d=>zrlx(:,:,jk)
92         type_array(jk)='T'
93         psgn_array(jk)=1.
94      END DO
95
96      !
97      IF( linit ) THEN              ! Metric coefficient (compute at the first call and saved in efact)
98         ALLOCATE( efact(jpi,jpj) , STAT=ierr )
99         IF( lk_mpp    )   CALL mpp_sum( ierr )
100         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'lim_hdf : unable to allocate arrays' )
101         DO jj = 2, jpjm1
102            DO ji = fs_2 , fs_jpim1   ! vector opt.
103               efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) * r1_e12t(ji,jj)
104            END DO
105         END DO
106         linit = .FALSE.
107      ENDIF
108      !                             ! Time integration parameters
109      !
110      zflu (jpi,: ) = 0._wp
111      zflv (jpi,: ) = 0._wp
112
113      DO jk=1 , isize
114         ztab0(:, : , jk ) = ptab(:,:,jk)      ! Arrays initialization
115         zdiv0(:, 1 , jk ) = 0._wp
116         zdiv0(:,jpj, jk ) = 0._wp
117         zdiv0(1,  :, jk ) = 0._wp
118         zdiv0(jpi,:, jk ) = 0._wp
119      END DO
120
121      zconv = 1._wp           !==  horizontal diffusion using a Crant-Nicholson scheme  ==!
122      iter  = 0
123      !
124      DO WHILE( MAXVAL(zconv(:)) > ( 2._wp * 1.e-04 ) .AND. iter <= its )   ! Sub-time step loop
125         !
126         iter = iter + 1                                 ! incrementation of the sub-time step number
127         !
128         DO jk = 1 , isize
129            jl = (jk-1) /( ihdf_vars+nlay_i)+1
130            IF (zconv(jk) > ( 2._wp * 1.e-04 )) THEN
131               DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
132                  DO ji = 1 , fs_jpim1   ! vector opt.
133                     zflu(ji,jj) = pahu3D(ji,jj,jl) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj,jk) - ptab(ji,jj,jk) )
134                     zflv(ji,jj) = pahv3D(ji,jj,jl) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1,jk) - ptab(ji,jj,jk) )
135                  END DO
136               END DO
137               !
138               DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes
139                  DO ji = fs_2 , fs_jpim1   ! vector opt.
140                     zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj)
141                  END DO
142               END DO
143               !
144               IF( iter == 1 )   zdiv0(:,:,jk) = zdiv(:,:)        ! save the 1st evaluation of the diffusive trend in zdiv0
145               !
146               DO jj = 2, jpjm1                                ! iterative evaluation
147                  DO ji = fs_2 , fs_jpim1   ! vector opt.
148                     zrlxint = (   ztab0(ji,jj,jk)    &
149                        &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj,jk) )   &
150                        &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj,jk) )                               &
151                        &      ) / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) )
152                     zrlx(ji,jj,jk) = ptab(ji,jj,jk) + zrelax * ( zrlxint - ptab(ji,jj,jk) )
153                  END DO
154               END DO
155            END IF
156
157         END DO
158
159         CALL lbc_lnk_multi( zrlx_array, type_array , psgn_array , isize ) ! Multiple interchange of all the variables
160         !
161         
162         IF ( MOD( iter-1 , nn_convfrq ) == 0 )  THEN   !Convergence test every nn_convfrq iterations (perf. optimization )
163            DO jk=1,isize
164               zconv(jk) = 0._wp                                   ! convergence test
165               DO jj = 2, jpjm1
166                  DO ji = fs_2, fs_jpim1
167                     zconv(jk) = MAX( zconv(jk), ABS( zrlx(ji,jj,jk) - ptab(ji,jj,jk) )  )
168                  END DO
169               END DO
170            END DO
171            IF( lk_mpp ) CALL mpp_max_multiple( zconv , isize )            ! max over the global domain for all the variables
172         ENDIF
173         !
174         DO jk=1,isize
175            ptab(:,:,jk) = zrlx(:,:,jk)
176         END DO
177         !
178      END DO                                       ! end of sub-time step loop
179
180     ! -----------------------
181      !!! final step (clem) !!!
182      DO jk = 1, isize
183         jl = (jk-1) /( ihdf_vars+nlay_i)+1
184         DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
185            DO ji = 1 , fs_jpim1   ! vector opt.
186               zflu(ji,jj) = pahu3D(ji,jj,jl) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj,jk) - ptab(ji,jj,jk) )
187               zflv(ji,jj) = pahv3D(ji,jj,jl) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1,jk) - ptab(ji,jj,jk) )
188            END DO
189         END DO
190         !
191         DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes
192            DO ji = fs_2 , fs_jpim1   ! vector opt.
193               zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj)
194               ptab(ji,jj,jk) = ztab0(ji,jj,jk) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj,jk) )
195            END DO
196         END DO
197      END DO
198
199      CALL lbc_lnk_multi( pt2d_array, type_array , psgn_array , isize ) ! Multiple interchange of all the variables
200
201      !!! final step (clem) !!!
202      ! -----------------------
203
204      IF(ln_ctl)   THEN
205         DO jk = 1 , isize
206            zrlx(:,:,jk) = ptab(:,:,jk) - ztab0(:,:,jk)
207            WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter
208            CALL prt_ctl( tab2d_1=zrlx(:,:,jk), clinfo1=charout )
209         END DO
210      ENDIF
211      !
212      CALL wrk_dealloc( jpi, jpj, isize, zrlx, zdiv0, ztab0 )
213      CALL wrk_dealloc( jpi, jpj, zflu, zflv, zdiv )
214
215      DEALLOCATE( zconv )
216      DEALLOCATE( pt2d_array , zrlx_array )
217      DEALLOCATE( type_array )
218      DEALLOCATE( psgn_array )
219      !
220   END SUBROUTINE lim_hdf
221
222
223   
224   SUBROUTINE lim_hdf_init
225      !!-------------------------------------------------------------------
226      !!                  ***  ROUTINE lim_hdf_init  ***
227      !!
228      !! ** Purpose : Initialisation of horizontal diffusion of sea-ice
229      !!
230      !! ** Method  : Read the namicehdf namelist
231      !!
232      !! ** input   : Namelist namicehdf
233      !!-------------------------------------------------------------------
234      INTEGER  ::   ios                 ! Local integer output status for namelist read
235      NAMELIST/namicehdf/ nn_convfrq 
236      !!-------------------------------------------------------------------
237      !
238      IF(lwp) THEN
239         WRITE(numout,*)
240         WRITE(numout,*) 'lim_hdf : Ice horizontal diffusion'
241         WRITE(numout,*) '~~~~~~~'
242      ENDIF
243      !
244      REWIND( numnam_ice_ref )              ! Namelist namicehdf in reference namelist : Ice horizontal diffusion
245      READ  ( numnam_ice_ref, namicehdf, IOSTAT = ios, ERR = 901)
246901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicehdf in reference namelist', lwp )
247
248      REWIND( numnam_ice_cfg )              ! Namelist namicehdf in configuration namelist : Ice horizontal diffusion
249      READ  ( numnam_ice_cfg, namicehdf, IOSTAT = ios, ERR = 902 )
250902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicehdf in configuration namelist', lwp )
251      IF(lwm) WRITE ( numoni, namicehdf )
252      !
253      IF(lwp) THEN                          ! control print
254         WRITE(numout,*)
255         WRITE(numout,*)'   Namelist of ice parameters for ice horizontal diffusion computation '
256         WRITE(numout,*)'      convergence check frequency of the Crant-Nicholson scheme   nn_convfrq   = ', nn_convfrq
257      ENDIF
258      !
259   END SUBROUTINE lim_hdf_init
260#else
261   !!----------------------------------------------------------------------
262   !!   Default option          Dummy module           NO LIM sea-ice model
263   !!----------------------------------------------------------------------
264#endif
265
266   !!======================================================================
267END MODULE limhdf
268
Note: See TracBrowser for help on using the repository browser.