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.
diapea.F90 in branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diapea.F90 @ 11038

Last change on this file since 11038 was 11038, checked in by rrenshaw, 5 years ago

code fixes
File size: 11.7 KB
RevLine 
[10958]1MODULE diapea
2   !!======================================================================
3   !!                       ***  MODULE  diapea  ***
4   !! Potential Energy Anomaly
5   !!======================================================================
6   !! History :  3.6  !  12/2016  (J Tinker)  Original code
7   !!----------------------------------------------------------------------
8   USE oce             ! ocean dynamics and tracers variables
9   USE dom_oce         ! ocean space and time domain
10   USE in_out_manager  ! I/O units
11   USE iom             ! I/0 library
12   USE wrk_nemo        ! working arrays
13   USE eosbn2          ! Equation of state - in situ and potential density
14   USE phycst          ! physical constant
15
16   IMPLICIT NONE
17   PRIVATE
18   
19   PUBLIC   dia_pea_init            ! routine called by nemogcm.F90
20   PUBLIC   dia_pea                 ! routine called by diawri.F90
21   REAL(wp), PUBLIC, SAVE, ALLOCATABLE,   DIMENSION(:,:)  ::   pea,peat,peas   
22   REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:,:)  ::   wgt_co_mat   ! Weighting array for proportion of grid shallower than cut off depth
23   REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:) ::   t_zmean, s_zmean  !Depth mean temperature and salinity: 2d fields
24   REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   t_zmean_mat, s_zmean_mat  !Depth mean temperature and salinity: 3d fields
25   REAL(wp) ::   zcutoff
26   LOGICAL , PUBLIC ::   ln_pea  ! region mean calculation
27   !!----------------------------------------------------------------------
28   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
29   !! $Id$
30   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
31   !!----------------------------------------------------------------------
32CONTAINS
33
34   SUBROUTINE dia_pea_init
35      ! Local variables
36      INTEGER :: ji,jj,jk  ! Dummy loop indices
37      REAL(wp) :: sumz,tmpsumz
38      INTEGER  :: ierr                ! error integer for IOM_get
39      INTEGER ::   ios                  ! Local integer output status for namelist read
40     
41      zcutoff = 200.!200m
42     
43      NAMELIST/nam_pea/ ln_pea
44     
45     
46      !
47      REWIND ( numnam_ref )              ! Read Namelist nam_diatmb in referdiatmbence namelist : TMB diagnostics
48      READ   ( numnam_ref, nam_pea, IOSTAT=ios, ERR= 901 )
49901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_pea in reference namelist', lwp )
50
51      REWIND( numnam_cfg )              ! Namelist nam_diatmb in configuration namelist  TMB diagnostics
52      READ  ( numnam_cfg, nam_pea, IOSTAT = ios, ERR = 902 )
53902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_pea in configuration namelist', lwp )
54      IF(lwm) WRITE ( numond, nam_pea )
55
56      IF(lwp) THEN                   ! Control print
57          WRITE(numout,*)
58          WRITE(numout,*) 'dia_pea_init : Output potential energy anomaly Diagnostics'
59          WRITE(numout,*) '~~~~~~~~~~~~'
60          WRITE(numout,*) 'Namelist nam_pea : set pea output '
[11038]61          WRITE(numout,*) 'Switch for pea diagnostics (T) or not (F)  ln_pea  = ', ln_pea
[10958]62      ENDIF
63     
64     
65   ALLOCATE( pea(jpi,jpj),  STAT= ierr )
66   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate pea array' )
67   
68   ALLOCATE( peat(jpi,jpj),  STAT= ierr )
69   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate peat array' )
70   
71   ALLOCATE( peas(jpi,jpj),  STAT= ierr )
72   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate peas array' )   
73     
74   ALLOCATE( t_zmean_mat(jpi,jpj,jpk),  STAT= ierr )
75   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate t_zmean_mat array' )
76   
77   ALLOCATE( s_zmean_mat(jpi,jpj,jpk),  STAT= ierr )
78   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate s_zmean_mat array' )
79   
80   ALLOCATE( t_zmean(jpi,jpj),  STAT= ierr )
81   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate t_zmean array' )
82   
83   ALLOCATE( s_zmean(jpi,jpj),  STAT= ierr )
84   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate s_zmean array' )
85   
86   ALLOCATE( wgt_co_mat(jpi,jpj,jpk),  STAT= ierr )
87   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate wgt_co_mat array' )
88   
89   
90   pea(:,:) = 0.
91   peat(:,:) = 0.
92   peas(:,:) = 0.
93   
94   if ( ln_pea ) THEN
95         
96      ! create wgt_co_mat mat, with the proportion of the grid (gdept_0) below cut off (200m)
97     
98      DO jj = 1,jpj
99          DO ji = 1,jpi
100              IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
101                sumz = 0.
102                DO jk = 1,jpk
103                  IF ( tmask(ji,jj,jk) == 1.0_wp ) THEN
104                    tmpsumz = sumz + e3t_n(ji,jj,jk)
105                    IF (sumz .ge. zcutoff) THEN
106                      ! Already too deep
107                      wgt_co_mat(ji,jj,jk) = 0.
108                    ELSE IF (tmpsumz .le. zcutoff) THEN
109                      ! Too shallow
110                      wgt_co_mat(ji,jj,jk) = 1.
111                    ELSE
112                      !proprotion of grid box above cut off depth
113                      wgt_co_mat(ji,jj,jk) = (zcutoff-Sumz)/e3t_n(ji,jj,jk)
114                    END IF
115                    sumz = tmpsumz
116                  endif
117                END DO
118               
119              ELSE
120                !if land, set to 0.
121                DO jk = 1,jpk
122                  wgt_co_mat(ji,jj,jk) = 0.
123                END DO
124                 
125              ENDIF
126          END DO
127       END DO
128   ENDIF
129   
130
131   END SUBROUTINE dia_pea_init
132   
133
134   SUBROUTINE dia_pea(kt)
135   INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
136   
137   INTEGER :: ji,jj,jk,ierr  ! Dummy loop indices
138   REAL(wp) :: tmpdenom, tmpnum, maxz
139   !rau0
140   
141   REAL(wp), ALLOCATABLE,   DIMENSION(:,:,:,:) ::   ts_pea_mat,ts_pea_mat_TS_mean,ts_pea_mat_S_mean
142   REAL(wp), ALLOCATABLE,   DIMENSION(:,:,:) ::   tmp_pea_rho,tmp_pea_TS_mean_rho,tmp_pea_S_mean_rho
143   REAL(wp) ::   int_y_pea,int_y_pea_t
144   
145   
146     
147   ALLOCATE( ts_pea_mat(jpi,jpj,jpk,2),  STAT= ierr )
148   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate ts_pea_mat array' )   
149   ALLOCATE( ts_pea_mat_TS_mean(jpi,jpj,jpk,2),  STAT= ierr )
150   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate ts_pea_mat_TS_mean array' )   
151   ALLOCATE( ts_pea_mat_S_mean(jpi,jpj,jpk,2),  STAT= ierr )
152   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate ts_pea_mat_S_mean array' )
153   ALLOCATE( tmp_pea_rho(jpi,jpj,jpk),  STAT= ierr )
154   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate tmp_pea_rho array' )
155   ALLOCATE( tmp_pea_TS_mean_rho(jpi,jpj,jpk),  STAT= ierr )
156   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate tmp_pea_TS_mean_rho array' )
157   ALLOCATE( tmp_pea_S_mean_rho(jpi,jpj,jpk),  STAT= ierr )
158   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate tmp_pea_S_mean_rho array' )
159
160    pea(:,:)=0.0d0
161    peaT(:,:)=0.0d0
162    peaS(:,:)=0.0d0
163   
164   
165   
166    !calculate the depth mean temperature and salinity of the upper 200m. Save this into a 3d array. Set the value where tmask=0 to be tsn.
167   
168     DO jj = 1,jpj
169         DO ji = 1,jpi
170            IF ( tmask(ji,jj,1) == 1.0_wp ) THEN  ! if a sea point.
171           
172             !Depth mean temperature
173             tmpdenom = 0.
174             tmpnum = 0.
175             DO jk = 1,jpk
176                IF ( tmask(ji,jj,jk) == 1.0_wp ) THEN
177                  tmpnum = tmpnum + (wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk)*tsn(ji,jj,jk,jp_tem))
178                  tmpdenom = tmpdenom + (wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk))
179                endif
180             END DO
181             t_zmean(ji,jj) = tmpnum/tmpdenom
182             
183             !Depth mean salinity
184             tmpdenom = 0.
185             tmpnum = 0.
186             DO jk = 1,jpk
187                IF ( tmask(ji,jj,jk) == 1.0_wp ) THEN
188                  tmpnum     = tmpnum + (wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk)*tsn(ji,jj,jk,jp_sal))
189                  tmpdenom = tmpdenom + (wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk))
190                endif
191             END DO
192             s_zmean(ji,jj) = tmpnum/tmpdenom
193             
194             !save into a 3d grid
195             DO jk = 1,jpk
196               IF ( tmask(ji,jj,jk) == 1.0_wp ) THEN
197                 t_zmean_mat(ji,jj,jk) = t_zmean(ji,jj)
198                 s_zmean_mat(ji,jj,jk) = s_zmean(ji,jj)
199               else
200                 t_zmean_mat(ji,jj,jk) = tsn(ji,jj,jk,jp_tem)
201                 s_zmean_mat(ji,jj,jk) = tsn(ji,jj,jk,jp_sal)
202               endif
203             END DO
204            else
205              t_zmean(ji,jj) = tsn(ji,jj,1,jp_tem)
206              s_zmean(ji,jj) = tsn(ji,jj,1,jp_sal)
207              DO jk = 1,jpk
208                t_zmean_mat(ji,jj,jk) = tsn(ji,jj,jk,jp_tem)
209                s_zmean_mat(ji,jj,jk) = tsn(ji,jj,jk,jp_sal)
210              END DO
211            endif
212           
213         END DO
214     END DO
215     
216   !Calculate the density from the depth varying, and depth average temperature and salinity
217   !-----------------------------
218   !-----------------------------
219   
220   ts_pea_mat(:,:,:,:) = tsn(:,:,:,:)
221   
222   ts_pea_mat_TS_mean(:,:,:,1) = t_zmean_mat(:,:,:)
223   ts_pea_mat_TS_mean(:,:,:,2) = s_zmean_mat(:,:,:)
224   
225   ts_pea_mat_S_mean(:,:,:,1) = t_zmean_mat(:,:,:)
226   ts_pea_mat_S_mean(:,:,:,2) = tsn(:,:,:,jp_sal)
227   
228   CALL eos ( ts_pea_mat,         tmp_pea_rho,         gdept_n(:,:,:) )
229   CALL eos ( ts_pea_mat_TS_mean, tmp_pea_TS_mean_rho, gdept_n(:,:,:) )
230   CALL eos ( ts_pea_mat_S_mean,  tmp_pea_S_mean_rho,  gdept_n(:,:,:) )
231   tmp_pea_rho = (tmp_pea_rho * rau0) + rau0
232   tmp_pea_TS_mean_rho = (tmp_pea_TS_mean_rho * rau0) + rau0
233   tmp_pea_S_mean_rho = (tmp_pea_S_mean_rho * rau0) + rau0
234   
235   
236   ! to test the density calculation
237   !CALL iom_put( "tmp_pea_rho" , tmp_pea_rho )                 ! pea
238   !CALL iom_put( "tmp_pea_TS_mean_rho" , tmp_pea_TS_mean_rho )                 ! pea
239   
240   
241   ! Caluclation of the PEA.
242    DO jj = 1,jpj
243      DO ji = 1,jpi
244        pea(ji,jj) = 0.
245        peat(ji,jj) = 0.
246        peas(ji,jj) = 0.
247        maxz = 0.
248        int_y_pea = 0.
249        int_y_pea_t = 0.
250        IF ( tmask(ji,jj,1) == 1.0_wp ) THEN ! for sea points
251       
252          ! the depth integrated calculation is summed up over the depths, and then divided by the depth
253          DO jk = 1,jpk
254            !for each level...
255           
256            IF ( tmask(ji,jj,jk) == 1.0_wp ) THEN ! if above the sea bed...
257              int_y_pea =   -((tmp_pea_TS_mean_rho(ji,jj,jk)) - (tmp_pea_rho(ji,jj,jk)))*9.81*gdept_n(ji,jj,jk)*wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk)
258              int_y_pea_t = -((tmp_pea_S_mean_rho(ji,jj,jk))  - (tmp_pea_rho(ji,jj,jk)))*9.81*gdept_n(ji,jj,jk)*wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk)
259            else
260              int_y_pea = 0.
261              int_y_pea_t = 0.
262            endif
263           
264            ! check that the sum is not NaN.
265            if ( int_y_pea .ne.  int_y_pea    ) int_y_pea = 0.
266            if ( int_y_pea_t .ne. int_y_pea_t )  int_y_pea_t = 0.
267            !if ( (int_y_pea*int_y_pea    ) .gt. 1.0e6 ) int_y_pea = 0.
268            !if ( (int_y_pea_t*int_y_pea_t) .gt. 1.0e6 ) int_y_pea_t = 0.
269           
270            pea(ji,jj) =  pea(ji,jj) + int_y_pea
271            peat(ji,jj) =  peat(ji,jj) + int_y_pea_t
272            maxz = maxz + (e3t_n(ji,jj,jk)*wgt_co_mat(ji,jj,jk))
273          enddo
274           
275           
276          !divide by the depth
277          pea(ji,jj) = pea(ji,jj)/maxz
278          peat(ji,jj) = peat(ji,jj)/maxz
279          peas(ji,jj) = pea(ji,jj) - peat(ji,jj)
280           
281           
282          else
283            pea(ji,jj) = 0.
284            peat(ji,jj) = 0.
285            peas(ji,jj) = 0.
286          endif
287       enddo
288    enddo
289!   
290     CALL iom_put( "pea" , pea )                 ! pea
291     CALL iom_put( "peat" , peat )               ! pea
292     CALL iom_put( "peas" , peas )               ! pea
293           
294
295    DEALLOCATE(ts_pea_mat,ts_pea_mat_TS_mean,ts_pea_mat_S_mean,tmp_pea_rho,tmp_pea_TS_mean_rho,tmp_pea_S_mean_rho)
296   
297   END SUBROUTINE dia_pea
298
299END MODULE diapea
Note: See TracBrowser for help on using the repository browser.