1 | MODULE 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 | !!---------------------------------------------------------------------- |
---|
32 | CONTAINS |
---|
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 ) |
---|
49 | 901 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 ) |
---|
53 | 902 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 ' |
---|
61 | WRITE(numout,*) 'Switch for pea diagnostics (T) or not (F) ln_diaregmean = ', ln_pea |
---|
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 | |
---|
299 | END MODULE diapea |
---|