source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90 @ 4528

Last change on this file since 4528 was 4161, checked in by cetlod, 7 years ago

dev_LOCEAN_2013 : merge in the 3rd dev branch dev_r4028_CNRS_LIM3, see ticket #1169

  • Property svn:keywords set to Id
File size: 6.9 KB
Line 
1MODULE limcons
2   !!======================================================================
3   !!                   ***  MODULE  limcons  ***
4   !! LIM-3 Sea Ice :   conservation check
5   !!======================================================================
6   !! History :   -   ! Original code from William H. Lipscomb, LANL
7   !!            3.0  ! 2004-06  (M. Vancoppenolle)   Energy Conservation
8   !!            4.0  ! 2011-02  (G. Madec)  add mpp considerations
9   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                      LIM-3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!    lim_cons     :   checks whether energy, mass and salt are conserved
15   !!----------------------------------------------------------------------
16   USE par_ice        ! LIM-3 parameter
17   USE ice            ! LIM-3 variables
18   USE dom_ice        ! LIM-3 domain
19   USE dom_oce        ! ocean domain
20   USE in_out_manager ! I/O manager
21   USE lib_mpp        ! MPP library
22   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   lim_column_sum
28   PUBLIC   lim_column_sum_energy
29   PUBLIC   lim_cons_check
30
31   !!----------------------------------------------------------------------
32   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
33   !! $Id$
34   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38   SUBROUTINE lim_column_sum( ksum, pin, pout )
39      !!-------------------------------------------------------------------
40      !!               ***  ROUTINE lim_column_sum ***
41      !!
42      !! ** Purpose : Compute the sum of xin over nsum categories
43      !!
44      !! ** Method  : Arithmetics
45      !!
46      !! ** Action  : Gets xin(ji,jj,jl) and computes xout(ji,jj)
47      !!---------------------------------------------------------------------
48      INTEGER                   , INTENT(in   ) ::   ksum   ! number of categories/layers
49      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pin    ! input field
50      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   pout   ! output field
51      !
52      INTEGER ::   jl   ! dummy loop indices
53      !!---------------------------------------------------------------------
54      !
55      pout(:,:) = pin(:,:,1)
56      DO jl = 2, ksum
57         pout(:,:) = pout(:,:) + pin(:,:,jl)
58      END DO
59      !
60   END SUBROUTINE lim_column_sum
61
62
63   SUBROUTINE lim_column_sum_energy( ksum, klay, pin, pout)
64      !!-------------------------------------------------------------------
65      !!               ***  ROUTINE lim_column_sum_energy ***
66      !!
67      !! ** Purpose : Compute the sum of xin over nsum categories
68      !!              and nlay layers
69      !!
70      !! ** Method  : Arithmetics
71      !!---------------------------------------------------------------------
72      INTEGER                               , INTENT(in   ) ::   ksum   !: number of categories
73      INTEGER                               , INTENT(in   ) ::   klay   !: number of vertical layers
74      REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl), INTENT(in   ) ::   pin   !: input field
75      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(  out) ::   pout   !: output field
76      !
77      INTEGER ::   jk, jl   ! dummy loop indices
78      !!---------------------------------------------------------------------
79      !
80      pout(:,:) = 0._wp
81      DO jl = 1, ksum
82         DO jk = 2, klay 
83            pout(:,:) = pout(:,:) + pin(:,:,jk,jl)
84         END DO
85      END DO
86      !
87   END SUBROUTINE lim_column_sum_energy
88
89
90   SUBROUTINE lim_cons_check( px1, px2, pmax_err, cd_fieldid )
91      !!-------------------------------------------------------------------
92      !!               ***  ROUTINE lim_cons_check ***
93      !!
94      !! ** Purpose : Test the conservation of a certain variable
95      !!              For each physical grid cell, check that initial
96      !!              and final values
97      !!              of a conserved field are equal to within a small value.
98      !!
99      !! ** Method  :
100      !!---------------------------------------------------------------------
101      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   px1          !: initial field
102      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   px2          !: final field
103      REAL(wp)                , INTENT(in   ) ::   pmax_err     !: max allowed error
104      CHARACTER(len=15)       , INTENT(in   ) ::   cd_fieldid   !: field identifyer
105      !
106      INTEGER  ::   ji, jj          ! dummy loop indices
107      INTEGER  ::   inb_error       ! number of g.c where there is a cons. error
108      LOGICAL  ::   llconserv_err   ! = .true. if conservation check failed
109      REAL(wp) ::   zmean_error     ! mean error on error points
110      !!---------------------------------------------------------------------
111      !
112      IF(lwp) WRITE(numout,*) ' lim_cons_check '
113      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
114
115      llconserv_err = .FALSE.
116      inb_error     = 0
117      zmean_error   = 0._wp
118      IF( MAXVAL( px2(:,:) - px1(:,:) ) > pmax_err )   llconserv_err = .TRUE.
119
120      IF( llconserv_err ) THEN
121         DO jj = 1, jpj 
122            DO ji = 1, jpi
123               IF( ABS( px2(ji,jj) - px1(ji,jj) ) > pmax_err ) THEN
124                  inb_error   = inb_error + 1
125                  zmean_error = zmean_error + ABS( px2(ji,jj) - px1(ji,jj) )
126                  !
127                  IF(lwp) THEN
128                     WRITE (numout,*) ' ALERTE 99 '
129                     WRITE (numout,*) ' Conservation error: ', cd_fieldid
130                     WRITE (numout,*) ' Point             : ', ji, jj 
131                     WRITE (numout,*) ' lat, lon          : ', gphit(ji,jj), glamt(ji,jj)
132                     WRITE (numout,*) ' Initial value     : ', px1(ji,jj)
133                     WRITE (numout,*) ' Final value       : ', px2(ji,jj)
134                     WRITE (numout,*) ' Difference        : ', px2(ji,jj) - px1(ji,jj)
135                  ENDIF
136               ENDIF
137            END DO
138         END DO
139         !
140      ENDIF
141      IF(lk_mpp)   CALL mpp_sum( inb_error   )
142      IF(lk_mpp)   CALL mpp_sum( zmean_error )
143      !
144      IF( inb_error > 0 .AND. lwp ) THEN
145         zmean_error = zmean_error / REAL( inb_error, wp )
146         WRITE(numout,*) ' Conservation check for : ', cd_fieldid
147         WRITE(numout,*) ' Number of error points : ', inb_error
148         WRITE(numout,*) ' Mean error on these pts: ', zmean_error
149      ENDIF
150      !
151   END SUBROUTINE lim_cons_check
152
153#else
154   !!----------------------------------------------------------------------
155   !!   Default option         Empty module            NO LIM sea-ice model
156   !!----------------------------------------------------------------------
157#endif
158   !!======================================================================
159END MODULE limcons
Note: See TracBrowser for help on using the repository browser.