Changeset 12377 for NEMO/trunk/src/OCE/DIA/diahth.F90
- Timestamp:
- 2020-02-12T15:39:06+01:00 (4 years ago)
- Location:
- NEMO/trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEAD ext/AGRIF5 ^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL
-
- Property svn:externals
-
NEMO/trunk/src/OCE/DIA/diahth.F90
r12276 r12377 40 40 41 41 42 !! * Substitutions 43 # include "do_loop_substitute.h90" 42 44 !!---------------------------------------------------------------------- 43 45 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 61 63 62 64 63 SUBROUTINE dia_hth( kt )65 SUBROUTINE dia_hth( kt, Kmm ) 64 66 !!--------------------------------------------------------------------- 65 67 !! *** ROUTINE dia_hth *** … … 82 84 !!------------------------------------------------------------------- 83 85 INTEGER, INTENT( in ) :: kt ! ocean time-step index 86 INTEGER, INTENT( in ) :: Kmm ! ocean time level index 84 87 !! 85 88 INTEGER :: ji, jj, jk ! dummy loop arguments … … 126 129 zdepinv(:,:) = 0._wp 127 130 zmaxdzT(:,:) = 0._wp 128 DO jj = 1, jpj 129 DO ji = 1, jpi 130 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 131 hth (ji,jj) = zztmp 132 zabs2 (ji,jj) = zztmp 133 ztm2 (ji,jj) = zztmp 134 zrho10_3(ji,jj) = zztmp 135 zpycn (ji,jj) = zztmp 136 END DO 137 END DO 131 DO_2D_11_11 132 zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 133 hth (ji,jj) = zztmp 134 zabs2 (ji,jj) = zztmp 135 ztm2 (ji,jj) = zztmp 136 zrho10_3(ji,jj) = zztmp 137 zpycn (ji,jj) = zztmp 138 END_2D 138 139 IF( nla10 > 1 ) THEN 139 DO jj = 1, jpj 140 DO ji = 1, jpi 141 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 142 zrho0_3(ji,jj) = zztmp 143 zrho0_1(ji,jj) = zztmp 144 END DO 145 END DO 140 DO_2D_11_11 141 zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 142 zrho0_3(ji,jj) = zztmp 143 zrho0_1(ji,jj) = zztmp 144 END_2D 146 145 ENDIF 147 146 148 147 ! Preliminary computation 149 148 ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 150 DO jj = 1, jpj 151 DO ji = 1, jpi 152 IF( tmask(ji,jj,nla10) == 1. ) THEN 153 zu = 1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80 * tsn(ji,jj,nla10,jp_sal) & 154 & - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) & 155 & - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal) 156 zv = 5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00 * tsn(ji,jj,nla10,jp_sal) & 157 & - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) 158 zut = 11.25 - 0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01 * tsn(ji,jj,nla10,jp_sal) 159 zvt = 38.00 - 0.750 * tsn(ji,jj,nla10,jp_tem) 160 zw = (zu + 0.698*zv) * (zu + 0.698*zv) 161 zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 162 ELSE 163 zdelr(ji,jj) = 0._wp 164 ENDIF 165 END DO 166 END DO 149 DO_2D_11_11 150 IF( tmask(ji,jj,nla10) == 1. ) THEN 151 zu = 1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80 * ts(ji,jj,nla10,jp_sal,Kmm) & 152 & - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm) & 153 & - 0.0100 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_sal,Kmm) 154 zv = 5891.00 + 38.000 * ts(ji,jj,nla10,jp_tem,Kmm) + 3.00 * ts(ji,jj,nla10,jp_sal,Kmm) & 155 & - 0.3750 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm) 156 zut = 11.25 - 0.149 * ts(ji,jj,nla10,jp_tem,Kmm) - 0.01 * ts(ji,jj,nla10,jp_sal,Kmm) 157 zvt = 38.00 - 0.750 * ts(ji,jj,nla10,jp_tem,Kmm) 158 zw = (zu + 0.698*zv) * (zu + 0.698*zv) 159 zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 160 ELSE 161 zdelr(ji,jj) = 0._wp 162 ENDIF 163 END_2D 167 164 168 165 ! ------------------------------------------------------------- ! … … 172 169 ! MLD: rho = rho(1) + zrho1 ! 173 170 ! ------------------------------------------------------------- ! 174 DO jk = jpkm1, 2, -1 ! loop from bottom to 2 175 DO jj = 1, jpj 176 DO ji = 1, jpi 177 ! 178 zzdep = gdepw_n(ji,jj,jk) 179 zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) & 180 & / zzdep * tmask(ji,jj,jk) ! vertical gradient of temperature (dT/dz) 181 zzdep = zzdep * tmask(ji,jj,1) 182 183 IF( zztmp > zmaxdzT(ji,jj) ) THEN 184 zmaxdzT(ji,jj) = zztmp 185 hth (ji,jj) = zzdep ! max and depth of dT/dz 186 ENDIF 187 188 IF( nla10 > 1 ) THEN 189 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1) ! delta rho(1) 190 IF( zztmp > zrho3 ) zrho0_3(ji,jj) = zzdep ! > 0.03 191 IF( zztmp > zrho1 ) zrho0_1(ji,jj) = zzdep ! > 0.01 192 ENDIF 193 END DO 194 END DO 195 END DO 171 DO_3DS_11_11( jpkm1, 2, -1 ) 172 ! 173 zzdep = gdepw(ji,jj,jk,Kmm) 174 zztmp = ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ) & 175 & / zzdep * tmask(ji,jj,jk) ! vertical gradient of temperature (dT/dz) 176 zzdep = zzdep * tmask(ji,jj,1) 177 178 IF( zztmp > zmaxdzT(ji,jj) ) THEN 179 zmaxdzT(ji,jj) = zztmp 180 hth (ji,jj) = zzdep ! max and depth of dT/dz 181 ENDIF 182 183 IF( nla10 > 1 ) THEN 184 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1) ! delta rho(1) 185 IF( zztmp > zrho3 ) zrho0_3(ji,jj) = zzdep ! > 0.03 186 IF( zztmp > zrho1 ) zrho0_1(ji,jj) = zzdep ! > 0.01 187 ENDIF 188 END_3D 196 189 197 190 CALL iom_put( 'mlddzt', hth ) ! depth of the thermocline … … 213 206 ! depth of temperature inversion ! 214 207 ! ------------------------------------------------------------- ! 215 DO jk = jpkm1, nlb10, -1 ! loop from bottom to nlb10 216 DO jj = 1, jpj 217 DO ji = 1, jpi 218 ! 219 zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1) 220 ! 221 zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem) ! - delta T(10m) 222 IF( ABS(zztmp) > ztem2 ) zabs2 (ji,jj) = zzdep ! abs > 0.2 223 IF( zztmp > ztem2 ) ztm2 (ji,jj) = zzdep ! > 0.2 224 zztmp = -zztmp ! delta T(10m) 225 IF( zztmp > ztinv(ji,jj) ) THEN ! temperature inversion 226 ztinv(ji,jj) = zztmp 227 zdepinv (ji,jj) = zzdep ! max value and depth 228 ENDIF 229 230 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10) ! delta rho(10m) 231 IF( zztmp > zrho3 ) zrho10_3(ji,jj) = zzdep ! > 0.03 232 IF( zztmp > zdelr(ji,jj) ) zpycn (ji,jj) = zzdep ! > equi. delta T(10m) - 0.2 233 ! 234 END DO 235 END DO 236 END DO 208 DO_3DS_11_11( jpkm1, nlb10, -1 ) 209 ! 210 zzdep = gdepw(ji,jj,jk,Kmm) * tmask(ji,jj,1) 211 ! 212 zztmp = ts(ji,jj,nla10,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ! - delta T(10m) 213 IF( ABS(zztmp) > ztem2 ) zabs2 (ji,jj) = zzdep ! abs > 0.2 214 IF( zztmp > ztem2 ) ztm2 (ji,jj) = zzdep ! > 0.2 215 zztmp = -zztmp ! delta T(10m) 216 IF( zztmp > ztinv(ji,jj) ) THEN ! temperature inversion 217 ztinv(ji,jj) = zztmp 218 zdepinv (ji,jj) = zzdep ! max value and depth 219 ENDIF 220 221 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10) ! delta rho(10m) 222 IF( zztmp > zrho3 ) zrho10_3(ji,jj) = zzdep ! > 0.03 223 IF( zztmp > zdelr(ji,jj) ) zpycn (ji,jj) = zzdep ! > equi. delta T(10m) - 0.2 224 ! 225 END_3D 237 226 238 227 CALL iom_put( 'mld_dt02', zabs2 ) ! MLD abs(delta t) - 0.2 … … 250 239 IF( iom_use ('20d') ) THEN ! depth of the 20 isotherm 251 240 ztem2 = 20. 252 CALL dia_hth_dep( ztem2, hd20 )241 CALL dia_hth_dep( Kmm, ztem2, hd20 ) 253 242 CALL iom_put( '20d', hd20 ) 254 243 ENDIF … … 256 245 IF( iom_use ('26d') ) THEN ! depth of the 26 isotherm 257 246 ztem2 = 26. 258 CALL dia_hth_dep( ztem2, hd26 )247 CALL dia_hth_dep( Kmm, ztem2, hd26 ) 259 248 CALL iom_put( '26d', hd26 ) 260 249 ENDIF … … 262 251 IF( iom_use ('28d') ) THEN ! depth of the 28 isotherm 263 252 ztem2 = 28. 264 CALL dia_hth_dep( ztem2, hd28 )253 CALL dia_hth_dep( Kmm, ztem2, hd28 ) 265 254 CALL iom_put( '28d', hd28 ) 266 255 ENDIF … … 271 260 IF( iom_use ('hc300') ) THEN 272 261 zzdep = 300. 273 CALL dia_hth_htc( zzdep, tsn(:,:,:,jp_tem), htc3 )262 CALL dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc3 ) 274 263 CALL iom_put( 'hc300', rau0_rcp * htc3 ) ! vertically integrated heat content (J/m2) 275 264 ENDIF … … 280 269 IF( iom_use ('hc700') ) THEN 281 270 zzdep = 700. 282 CALL dia_hth_htc( zzdep, tsn(:,:,:,jp_tem), htc7 )271 CALL dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc7 ) 283 272 CALL iom_put( 'hc700', rau0_rcp * htc7 ) ! vertically integrated heat content (J/m2) 284 273 … … 290 279 IF( iom_use ('hc2000') ) THEN 291 280 zzdep = 2000. 292 CALL dia_hth_htc( zzdep, tsn(:,:,:,jp_tem), htc20 )281 CALL dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc20 ) 293 282 CALL iom_put( 'hc2000', rau0_rcp * htc20 ) ! vertically integrated heat content (J/m2) 294 283 ENDIF … … 301 290 END SUBROUTINE dia_hth 302 291 303 SUBROUTINE dia_hth_dep( ptem, pdept ) 304 ! 292 SUBROUTINE dia_hth_dep( Kmm, ptem, pdept ) 293 ! 294 INTEGER , INTENT(in) :: Kmm ! ocean time level index 305 295 REAL(wp), INTENT(in) :: ptem 306 296 REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdept … … 314 304 ! --------------------------------------- ! 315 305 iktem(:,:) = 1 316 DO jk = 1, jpkm1 ! beware temperature is not always decreasing with depth => loop from top to bottom 317 DO jj = 1, jpj 318 DO ji = 1, jpi 319 zztmp = tsn(ji,jj,jk,jp_tem) 320 IF( zztmp >= ptem ) iktem(ji,jj) = jk 321 END DO 322 END DO 323 END DO 306 DO_3D_11_11( 1, jpkm1 ) 307 zztmp = ts(ji,jj,jk,jp_tem,Kmm) 308 IF( zztmp >= ptem ) iktem(ji,jj) = jk 309 END_3D 324 310 325 311 ! ------------------------------- ! 326 312 ! Depth of ptem isotherm ! 327 313 ! ------------------------------- ! 328 DO jj = 1, jpj 329 DO ji = 1, jpi 330 ! 331 zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean bottom 332 ! 333 iid = iktem(ji,jj) 334 IF( iid /= 1 ) THEN 335 zztmp = gdept_n(ji,jj,iid ) & ! linear interpolation 336 & + ( gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid) ) & 337 & * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem) ) & 338 & / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 339 pdept(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1) ! bound by the ocean depth 340 ELSE 341 pdept(ji,jj) = 0._wp 342 ENDIF 343 END DO 344 END DO 314 DO_2D_11_11 315 ! 316 zzdep = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ! depth of the ocean bottom 317 ! 318 iid = iktem(ji,jj) 319 IF( iid /= 1 ) THEN 320 zztmp = gdept(ji,jj,iid ,Kmm) & ! linear interpolation 321 & + ( gdept(ji,jj,iid+1,Kmm) - gdept(ji,jj,iid,Kmm) ) & 322 & * ( 20.*tmask(ji,jj,iid+1) - ts(ji,jj,iid,jp_tem,Kmm) ) & 323 & / ( ts(ji,jj,iid+1,jp_tem,Kmm) - ts(ji,jj,iid,jp_tem,Kmm) + (1.-tmask(ji,jj,1)) ) 324 pdept(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1) ! bound by the ocean depth 325 ELSE 326 pdept(ji,jj) = 0._wp 327 ENDIF 328 END_2D 345 329 ! 346 330 END SUBROUTINE dia_hth_dep 347 331 348 332 349 SUBROUTINE dia_hth_htc( pdep, ptn, phtc ) 350 ! 333 SUBROUTINE dia_hth_htc( Kmm, pdep, pt, phtc ) 334 ! 335 INTEGER , INTENT(in) :: Kmm ! ocean time level index 351 336 REAL(wp), INTENT(in) :: pdep ! depth over the heat content 352 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pt n337 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pt 353 338 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phtc 354 339 ! … … 361 346 362 347 IF( .NOT. ln_linssh ) THEN ; zthick(:,:) = 0._wp ; phtc(:,:) = 0._wp 363 ELSE ; zthick(:,:) = ssh n(:,:) ; phtc(:,:) = ptn(:,:,1) * sshn(:,:) * tmask(:,:,1)348 ELSE ; zthick(:,:) = ssh(:,:,Kmm) ; phtc(:,:) = pt(:,:,1) * ssh(:,:,Kmm) * tmask(:,:,1) 364 349 ENDIF 365 350 ! 366 351 ilevel(:,:) = 1 367 DO jk = 2, jpkm1 368 DO jj = 1, jpj 369 DO ji = 1, jpi 370 IF( ( gdept_n(ji,jj,jk) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN 371 ilevel(ji,jj) = jk 372 zthick(ji,jj) = zthick(ji,jj) + e3t_n(ji,jj,jk) 373 phtc (ji,jj) = phtc (ji,jj) + e3t_n(ji,jj,jk) * ptn(ji,jj,jk) 374 ENDIF 375 ENDDO 376 ENDDO 377 ENDDO 378 ! 379 DO jj = 1, jpj 380 DO ji = 1, jpi 381 ik = ilevel(ji,jj) 382 zthick(ji,jj) = pdep - zthick(ji,jj) ! remaining thickness to reach depht pdep 383 phtc(ji,jj) = phtc(ji,jj) + ptn(ji,jj,ik+1) * MIN( e3t_n(ji,jj,ik+1), zthick(ji,jj) ) & 384 * tmask(ji,jj,ik+1) 385 END DO 386 ENDDO 352 DO_3D_11_11( 2, jpkm1 ) 353 IF( ( gdept(ji,jj,jk,Kmm) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN 354 ilevel(ji,jj) = jk 355 zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm) 356 phtc (ji,jj) = phtc (ji,jj) + e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk) 357 ENDIF 358 END_3D 359 ! 360 DO_2D_11_11 361 ik = ilevel(ji,jj) 362 zthick(ji,jj) = pdep - zthick(ji,jj) ! remaining thickness to reach depht pdep 363 phtc(ji,jj) = phtc(ji,jj) + pt(ji,jj,ik+1) * MIN( e3t(ji,jj,ik+1,Kmm), zthick(ji,jj) ) & 364 * tmask(ji,jj,ik+1) 365 END_2D 387 366 ! 388 367 !
Note: See TracChangeset
for help on using the changeset viewer.