- Timestamp:
- 2019-12-11T17:15:54+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DIA/diahth.F90
r11949 r12193 11 11 !! 3.2 ! 2009-07 (S. Masson) hc300 bugfix + cleaning + add new diag 12 12 !!---------------------------------------------------------------------- 13 #if defined key_diahth14 !!----------------------------------------------------------------------15 !! 'key_diahth' : thermocline depth diag.16 !!----------------------------------------------------------------------17 13 !! dia_hth : Compute varius diagnostics associated with the mixed layer 18 14 !!---------------------------------------------------------------------- … … 32 28 PUBLIC dia_hth_alloc ! routine called by nemogcm.F90 33 29 34 LOGICAL , PUBLIC, PARAMETER :: lk_diahth = .TRUE.!: thermocline-20d depths flag30 LOGICAL, SAVE :: l_hth !: thermocline-20d depths flag 35 31 36 32 ! note: following variables should move to local variables once iom_put is always used 37 33 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hth !: depth of the max vertical temperature gradient [m] 38 34 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hd20 !: depth of 20 C isotherm [m] 35 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hd26 !: depth of 26 C isotherm [m] 39 36 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hd28 !: depth of 28 C isotherm [m] 40 37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: htc3 !: heat content of first 300 m [W] 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: htc7 !: heat content of first 700 m [W] 39 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: htc20 !: heat content of first 2000 m [W] 40 41 41 42 42 !!---------------------------------------------------------------------- … … 52 52 !!--------------------------------------------------------------------- 53 53 ! 54 ALLOCATE( hth(jpi,jpj), hd20(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), STAT=dia_hth_alloc ) 54 ALLOCATE( hth(jpi,jpj), hd20(jpi,jpj), hd26(jpi,jpj), hd28(jpi,jpj), & 55 & htc3(jpi,jpj), htc7(jpi,jpj), htc20(jpi,jpj), STAT=dia_hth_alloc ) 55 56 ! 56 57 CALL mpp_sum ( 'diahth', dia_hth_alloc ) … … 83 84 INTEGER, INTENT( in ) :: Kmm ! ocean time level index 84 85 !! 85 INTEGER :: ji, jj, jk ! dummy loop arguments 86 INTEGER :: iid, ilevel ! temporary integers 87 INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ik20, ik28 ! levels 88 REAL(wp) :: zavt5 = 5.e-4_wp ! Kz criterion for the turbocline depth 89 REAL(wp) :: zrho3 = 0.03_wp ! density criterion for mixed layer depth 90 REAL(wp) :: zrho1 = 0.01_wp ! density criterion for mixed layer depth 91 REAL(wp) :: ztem2 = 0.2_wp ! temperature criterion for mixed layer depth 92 REAL(wp) :: zthick_0, zcoef ! temporary scalars 93 REAL(wp) :: zztmp, zzdep ! temporary scalars inside do loop 94 REAL(wp) :: zu, zv, zw, zut, zvt ! temporary workspace 95 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zabs2 ! MLD: abs( tn - tn(10m) ) = ztem2 96 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ztm2 ! Top of thermocline: tn = tn(10m) - ztem2 97 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zrho10_3 ! MLD: rho = rho10m + zrho3 98 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zpycn ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 99 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ztinv ! max of temperature inversion 100 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdepinv ! depth of temperature inversion 101 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zrho0_3 ! MLD rho = rho(surf) = 0.03 102 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zrho0_1 ! MLD rho = rho(surf) = 0.01 103 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zmaxdzT ! max of dT/dz 104 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zthick ! vertical integration thickness 105 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdelr ! delta rho equivalent to deltaT = 0.2 86 INTEGER :: ji, jj, jk ! dummy loop arguments 87 REAL(wp) :: zrho3 = 0.03_wp ! density criterion for mixed layer depth 88 REAL(wp) :: zrho1 = 0.01_wp ! density criterion for mixed layer depth 89 REAL(wp) :: ztem2 = 0.2_wp ! temperature criterion for mixed layer depth 90 REAL(wp) :: zztmp, zzdep ! temporary scalars inside do loop 91 REAL(wp) :: zu, zv, zw, zut, zvt ! temporary workspace 92 REAL(wp), DIMENSION(jpi,jpj) :: zabs2 ! MLD: abs( tn - tn(10m) ) = ztem2 93 REAL(wp), DIMENSION(jpi,jpj) :: ztm2 ! Top of thermocline: tn = tn(10m) - ztem2 94 REAL(wp), DIMENSION(jpi,jpj) :: zrho10_3 ! MLD: rho = rho10m + zrho3 95 REAL(wp), DIMENSION(jpi,jpj) :: zpycn ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 96 REAL(wp), DIMENSION(jpi,jpj) :: ztinv ! max of temperature inversion 97 REAL(wp), DIMENSION(jpi,jpj) :: zdepinv ! depth of temperature inversion 98 REAL(wp), DIMENSION(jpi,jpj) :: zrho0_3 ! MLD rho = rho(surf) = 0.03 99 REAL(wp), DIMENSION(jpi,jpj) :: zrho0_1 ! MLD rho = rho(surf) = 0.01 100 REAL(wp), DIMENSION(jpi,jpj) :: zmaxdzT ! max of dT/dz 101 REAL(wp), DIMENSION(jpi,jpj) :: zdelr ! delta rho equivalent to deltaT = 0.2 106 102 !!---------------------------------------------------------------------- 107 103 IF( ln_timing ) CALL timing_start('dia_hth') 108 104 109 105 IF( kt == nit000 ) THEN 106 l_hth = .FALSE. 107 IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) .OR. & 108 & iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) .OR. iom_use( 'mldr10_3' ) .OR. & 109 & iom_use( '20d' ) .OR. iom_use( '26d' ) .OR. iom_use( '28d' ) .OR. & 110 & iom_use( 'hc300' ) .OR. iom_use( 'hc700' ) .OR. iom_use( 'hc2000' ) .OR. & 111 & iom_use( 'pycndep' ) .OR. iom_use( 'tinv' ) .OR. iom_use( 'depti' ) ) l_hth = .TRUE. 110 112 ! ! allocate dia_hth array 111 IF( dia_hth_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard arrays' ) 112 113 IF(.NOT. ALLOCATED(ik20) ) THEN 114 ALLOCATE(ik20(jpi,jpj), ik28(jpi,jpj), & 115 & zabs2(jpi,jpj), & 116 & ztm2(jpi,jpj), & 117 & zrho10_3(jpi,jpj),& 118 & zpycn(jpi,jpj), & 119 & ztinv(jpi,jpj), & 120 & zdepinv(jpi,jpj), & 121 & zrho0_3(jpi,jpj), & 122 & zrho0_1(jpi,jpj), & 123 & zmaxdzT(jpi,jpj), & 124 & zthick(jpi,jpj), & 125 & zdelr(jpi,jpj), STAT=ji) 126 CALL mpp_sum('diahth', ji) 127 IF( ji /= 0 ) CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard ocean arrays' ) 128 END IF 129 130 IF(lwp) WRITE(numout,*) 131 IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth' 132 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 133 IF(lwp) WRITE(numout,*) 113 IF( l_hth ) THEN 114 IF( dia_hth_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard arrays' ) 115 IF(lwp) WRITE(numout,*) 116 IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth' 117 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 118 IF(lwp) WRITE(numout,*) 119 ENDIF 134 120 ENDIF 135 121 136 ! initialization 137 ztinv (:,:) = 0._wp 138 zdepinv(:,:) = 0._wp 139 zmaxdzT(:,:) = 0._wp 140 DO jj = 1, jpj 141 DO ji = 1, jpi 142 zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 143 hth (ji,jj) = zztmp 144 zabs2 (ji,jj) = zztmp 145 ztm2 (ji,jj) = zztmp 146 zrho10_3(ji,jj) = zztmp 147 zpycn (ji,jj) = zztmp 148 END DO 149 END DO 150 IF( nla10 > 1 ) THEN 151 DO jj = 1, jpj 152 DO ji = 1, jpi 153 zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 154 zrho0_3(ji,jj) = zztmp 155 zrho0_1(ji,jj) = zztmp 156 END DO 157 END DO 122 IF( l_hth ) THEN 123 ! 124 IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) ) THEN 125 ! initialization 126 ztinv (:,:) = 0._wp 127 zdepinv(:,:) = 0._wp 128 zmaxdzT(:,:) = 0._wp 129 DO jj = 1, jpj 130 DO ji = 1, jpi 131 zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 132 hth (ji,jj) = zztmp 133 zabs2 (ji,jj) = zztmp 134 ztm2 (ji,jj) = zztmp 135 zrho10_3(ji,jj) = zztmp 136 zpycn (ji,jj) = zztmp 137 END DO 138 END DO 139 IF( nla10 > 1 ) THEN 140 DO jj = 1, jpj 141 DO ji = 1, jpi 142 zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 143 zrho0_3(ji,jj) = zztmp 144 zrho0_1(ji,jj) = zztmp 145 END DO 146 END DO 147 ENDIF 148 149 ! Preliminary computation 150 ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 151 DO jj = 1, jpj 152 DO ji = 1, jpi 153 IF( tmask(ji,jj,nla10) == 1. ) THEN 154 zu = 1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80 * ts(ji,jj,nla10,jp_sal,Kmm) & 155 & - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm) & 156 & - 0.0100 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_sal,Kmm) 157 zv = 5891.00 + 38.000 * ts(ji,jj,nla10,jp_tem,Kmm) + 3.00 * ts(ji,jj,nla10,jp_sal,Kmm) & 158 & - 0.3750 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm) 159 zut = 11.25 - 0.149 * ts(ji,jj,nla10,jp_tem,Kmm) - 0.01 * ts(ji,jj,nla10,jp_sal,Kmm) 160 zvt = 38.00 - 0.750 * ts(ji,jj,nla10,jp_tem,Kmm) 161 zw = (zu + 0.698*zv) * (zu + 0.698*zv) 162 zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 163 ELSE 164 zdelr(ji,jj) = 0._wp 165 ENDIF 166 END DO 167 END DO 168 169 ! ------------------------------------------------------------- ! 170 ! thermocline depth: strongest vertical gradient of temperature ! 171 ! turbocline depth (mixing layer depth): avt = zavt5 ! 172 ! MLD: rho = rho(1) + zrho3 ! 173 ! MLD: rho = rho(1) + zrho1 ! 174 ! ------------------------------------------------------------- ! 175 DO jk = jpkm1, 2, -1 ! loop from bottom to 2 176 DO jj = 1, jpj 177 DO ji = 1, jpi 178 ! 179 zzdep = gdepw(ji,jj,jk,Kmm) 180 zztmp = ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ) & 181 & / zzdep * tmask(ji,jj,jk) ! vertical gradient of temperature (dT/dz) 182 zzdep = zzdep * tmask(ji,jj,1) 183 184 IF( zztmp > zmaxdzT(ji,jj) ) THEN 185 zmaxdzT(ji,jj) = zztmp 186 hth (ji,jj) = zzdep ! max and depth of dT/dz 187 ENDIF 188 189 IF( nla10 > 1 ) THEN 190 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1) ! delta rho(1) 191 IF( zztmp > zrho3 ) zrho0_3(ji,jj) = zzdep ! > 0.03 192 IF( zztmp > zrho1 ) zrho0_1(ji,jj) = zzdep ! > 0.01 193 ENDIF 194 END DO 195 END DO 196 END DO 197 198 CALL iom_put( 'mlddzt', hth ) ! depth of the thermocline 199 IF( nla10 > 1 ) THEN 200 CALL iom_put( 'mldr0_3', zrho0_3 ) ! MLD delta rho(surf) = 0.03 201 CALL iom_put( 'mldr0_1', zrho0_1 ) ! MLD delta rho(surf) = 0.01 202 ENDIF 203 ! 204 ENDIF 205 ! 206 IF( iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) .OR. iom_use( 'mldr10_3' ) .OR. & 207 & iom_use( 'pycndep' ) .OR. iom_use( 'tinv' ) .OR. iom_use( 'depti' ) ) THEN 208 ! ------------------------------------------------------------- ! 209 ! MLD: abs( tn - tn(10m) ) = ztem2 ! 210 ! Top of thermocline: tn = tn(10m) - ztem2 ! 211 ! MLD: rho = rho10m + zrho3 ! 212 ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) ! 213 ! temperature inversion: max( 0, max of tn - tn(10m) ) ! 214 ! depth of temperature inversion ! 215 ! ------------------------------------------------------------- ! 216 DO jk = jpkm1, nlb10, -1 ! loop from bottom to nlb10 217 DO jj = 1, jpj 218 DO ji = 1, jpi 219 ! 220 zzdep = gdepw(ji,jj,jk,Kmm) * tmask(ji,jj,1) 221 ! 222 zztmp = ts(ji,jj,nla10,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ! - delta T(10m) 223 IF( ABS(zztmp) > ztem2 ) zabs2 (ji,jj) = zzdep ! abs > 0.2 224 IF( zztmp > ztem2 ) ztm2 (ji,jj) = zzdep ! > 0.2 225 zztmp = -zztmp ! delta T(10m) 226 IF( zztmp > ztinv(ji,jj) ) THEN ! temperature inversion 227 ztinv(ji,jj) = zztmp 228 zdepinv (ji,jj) = zzdep ! max value and depth 229 ENDIF 230 231 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10) ! delta rho(10m) 232 IF( zztmp > zrho3 ) zrho10_3(ji,jj) = zzdep ! > 0.03 233 IF( zztmp > zdelr(ji,jj) ) zpycn (ji,jj) = zzdep ! > equi. delta T(10m) - 0.2 234 ! 235 END DO 236 END DO 237 END DO 238 239 CALL iom_put( 'mld_dt02', zabs2 ) ! MLD abs(delta t) - 0.2 240 CALL iom_put( 'topthdep', ztm2 ) ! T(10) - 0.2 241 CALL iom_put( 'mldr10_3', zrho10_3 ) ! MLD delta rho(10m) = 0.03 242 CALL iom_put( 'pycndep' , zpycn ) ! MLD delta rho equi. delta T(10m) = 0.2 243 CALL iom_put( 'tinv' , ztinv ) ! max. temp. inv. (t10 ref) 244 CALL iom_put( 'depti' , zdepinv ) ! depth of max. temp. inv. (t10 ref) 245 ! 246 ENDIF 247 248 ! ------------------------------- ! 249 ! Depth of 20C/26C/28C isotherm ! 250 ! ------------------------------- ! 251 IF( iom_use ('20d') ) THEN ! depth of the 20 isotherm 252 ztem2 = 20. 253 CALL dia_hth_dep( Kmm, ztem2, hd20 ) 254 CALL iom_put( '20d', hd20 ) 255 ENDIF 256 ! 257 IF( iom_use ('26d') ) THEN ! depth of the 26 isotherm 258 ztem2 = 26. 259 CALL dia_hth_dep( Kmm, ztem2, hd26 ) 260 CALL iom_put( '26d', hd26 ) 261 ENDIF 262 ! 263 IF( iom_use ('28d') ) THEN ! depth of the 28 isotherm 264 ztem2 = 28. 265 CALL dia_hth_dep( Kmm, ztem2, hd28 ) 266 CALL iom_put( '28d', hd28 ) 267 ENDIF 268 269 ! ----------------------------- ! 270 ! Heat content of first 300 m ! 271 ! ----------------------------- ! 272 IF( iom_use ('hc300') ) THEN 273 zzdep = 300. 274 CALL dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc3 ) 275 CALL iom_put( 'hc300', rau0_rcp * htc3 ) ! vertically integrated heat content (J/m2) 276 ENDIF 277 ! 278 ! ----------------------------- ! 279 ! Heat content of first 700 m ! 280 ! ----------------------------- ! 281 IF( iom_use ('hc700') ) THEN 282 zzdep = 700. 283 CALL dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc7 ) 284 CALL iom_put( 'hc700', rau0_rcp * htc7 ) ! vertically integrated heat content (J/m2) 285 286 ENDIF 287 ! 288 ! ----------------------------- ! 289 ! Heat content of first 2000 m ! 290 ! ----------------------------- ! 291 IF( iom_use ('hc2000') ) THEN 292 zzdep = 2000. 293 CALL dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc20 ) 294 CALL iom_put( 'hc2000', rau0_rcp * htc20 ) ! vertically integrated heat content (J/m2) 295 ENDIF 296 ! 158 297 ENDIF 298 299 ! 300 IF( ln_timing ) CALL timing_stop('dia_hth') 301 ! 302 END SUBROUTINE dia_hth 303 304 SUBROUTINE dia_hth_dep( Kmm, ptem, pdept ) 305 ! 306 INTEGER , INTENT(in) :: Kmm ! ocean time level index 307 REAL(wp), INTENT(in) :: ptem 308 REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdept 309 ! 310 INTEGER :: ji, jj, jk, iid 311 REAL(wp) :: zztmp, zzdep 312 INTEGER, DIMENSION(jpi,jpj) :: iktem 159 313 160 ! Preliminary computation 161 ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 162 DO jj = 1, jpj 163 DO ji = 1, jpi 164 IF( tmask(ji,jj,nla10) == 1. ) THEN 165 zu = 1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80 * ts(ji,jj,nla10,jp_sal,Kmm) & 166 & - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm) & 167 & - 0.0100 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_sal,Kmm) 168 zv = 5891.00 + 38.000 * ts(ji,jj,nla10,jp_tem,Kmm) + 3.00 * ts(ji,jj,nla10,jp_sal,Kmm) & 169 & - 0.3750 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm) 170 zut = 11.25 - 0.149 * ts(ji,jj,nla10,jp_tem,Kmm) - 0.01 * ts(ji,jj,nla10,jp_sal,Kmm) 171 zvt = 38.00 - 0.750 * ts(ji,jj,nla10,jp_tem,Kmm) 172 zw = (zu + 0.698*zv) * (zu + 0.698*zv) 173 zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 174 ELSE 175 zdelr(ji,jj) = 0._wp 176 ENDIF 177 END DO 178 END DO 179 180 ! ------------------------------------------------------------- ! 181 ! thermocline depth: strongest vertical gradient of temperature ! 182 ! turbocline depth (mixing layer depth): avt = zavt5 ! 183 ! MLD: rho = rho(1) + zrho3 ! 184 ! MLD: rho = rho(1) + zrho1 ! 185 ! ------------------------------------------------------------- ! 186 DO jk = jpkm1, 2, -1 ! loop from bottom to 2 187 DO jj = 1, jpj 188 DO ji = 1, jpi 189 ! 190 zzdep = gdepw(ji,jj,jk,Kmm) 191 zztmp = ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ) / zzdep * tmask(ji,jj,jk) ! vertical gradient of temperature (dT/dz) 192 zzdep = zzdep * tmask(ji,jj,1) 193 194 IF( zztmp > zmaxdzT(ji,jj) ) THEN 195 zmaxdzT(ji,jj) = zztmp ; hth (ji,jj) = zzdep ! max and depth of dT/dz 196 ENDIF 197 198 IF( nla10 > 1 ) THEN 199 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1) ! delta rho(1) 200 IF( zztmp > zrho3 ) zrho0_3(ji,jj) = zzdep ! > 0.03 201 IF( zztmp > zrho1 ) zrho0_1(ji,jj) = zzdep ! > 0.01 202 ENDIF 203 204 END DO 205 END DO 206 END DO 207 208 CALL iom_put( "mlddzt", hth ) ! depth of the thermocline 209 IF( nla10 > 1 ) THEN 210 CALL iom_put( "mldr0_3", zrho0_3 ) ! MLD delta rho(surf) = 0.03 211 CALL iom_put( "mldr0_1", zrho0_1 ) ! MLD delta rho(surf) = 0.01 212 ENDIF 213 214 ! ------------------------------------------------------------- ! 215 ! MLD: abs( tn - tn(10m) ) = ztem2 ! 216 ! Top of thermocline: tn = tn(10m) - ztem2 ! 217 ! MLD: rho = rho10m + zrho3 ! 218 ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) ! 219 ! temperature inversion: max( 0, max of tn - tn(10m) ) ! 220 ! depth of temperature inversion ! 221 ! ------------------------------------------------------------- ! 222 DO jk = jpkm1, nlb10, -1 ! loop from bottom to nlb10 223 DO jj = 1, jpj 224 DO ji = 1, jpi 225 ! 226 zzdep = gdepw(ji,jj,jk,Kmm) * tmask(ji,jj,1) 227 ! 228 zztmp = ts(ji,jj,nla10,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ! - delta T(10m) 229 IF( ABS(zztmp) > ztem2 ) zabs2 (ji,jj) = zzdep ! abs > 0.2 230 IF( zztmp > ztem2 ) ztm2 (ji,jj) = zzdep ! > 0.2 231 zztmp = -zztmp ! delta T(10m) 232 IF( zztmp > ztinv(ji,jj) ) THEN ! temperature inversion 233 ztinv(ji,jj) = zztmp ; zdepinv (ji,jj) = zzdep ! max value and depth 234 ENDIF 235 236 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10) ! delta rho(10m) 237 IF( zztmp > zrho3 ) zrho10_3(ji,jj) = zzdep ! > 0.03 238 IF( zztmp > zdelr(ji,jj) ) zpycn (ji,jj) = zzdep ! > equi. delta T(10m) - 0.2 239 ! 240 END DO 241 END DO 242 END DO 243 244 CALL iom_put( "mld_dt02", zabs2 ) ! MLD abs(delta t) - 0.2 245 CALL iom_put( "topthdep", ztm2 ) ! T(10) - 0.2 246 CALL iom_put( "mldr10_3", zrho10_3 ) ! MLD delta rho(10m) = 0.03 247 CALL iom_put( "pycndep" , zpycn ) ! MLD delta rho equi. delta T(10m) = 0.2 248 CALL iom_put( "tinv" , ztinv ) ! max. temp. inv. (t10 ref) 249 CALL iom_put( "depti" , zdepinv ) ! depth of max. temp. inv. (t10 ref) 250 251 252 ! ----------------------------------- ! 253 ! search deepest level above 20C/28C ! 254 ! ----------------------------------- ! 255 ik20(:,:) = 1 256 ik28(:,:) = 1 314 ! --------------------------------------- ! 315 ! search deepest level above ptem ! 316 ! --------------------------------------- ! 317 iktem(:,:) = 1 257 318 DO jk = 1, jpkm1 ! beware temperature is not always decreasing with depth => loop from top to bottom 258 319 DO jj = 1, jpj 259 320 DO ji = 1, jpi 260 321 zztmp = ts(ji,jj,jk,jp_tem,Kmm) 261 IF( zztmp >= 20. ) ik20(ji,jj) = jk 262 IF( zztmp >= 28. ) ik28(ji,jj) = jk 322 IF( zztmp >= ptem ) iktem(ji,jj) = jk 263 323 END DO 264 324 END DO 265 325 END DO 266 326 267 ! --------------------------- !268 ! Depth of 20C/28C isotherm!269 ! --------------------------- !327 ! ------------------------------- ! 328 ! Depth of ptem isotherm ! 329 ! ------------------------------- ! 270 330 DO jj = 1, jpj 271 331 DO ji = 1, jpi 272 332 ! 273 zzdep = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ! depth of the o ean bottom333 zzdep = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ! depth of the ocean bottom 274 334 ! 275 iid = ik 20(ji,jj)335 iid = iktem(ji,jj) 276 336 IF( iid /= 1 ) THEN 277 zztmp =gdept(ji,jj,iid ,Kmm) & ! linear interpolation337 zztmp = gdept(ji,jj,iid ,Kmm) & ! linear interpolation 278 338 & + ( gdept(ji,jj,iid+1,Kmm) - gdept(ji,jj,iid,Kmm) ) & 279 339 & * ( 20.*tmask(ji,jj,iid+1) - ts(ji,jj,iid,jp_tem,Kmm) ) & 280 340 & / ( ts(ji,jj,iid+1,jp_tem,Kmm) - ts(ji,jj,iid,jp_tem,Kmm) + (1.-tmask(ji,jj,1)) ) 281 hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1) ! bound by the ocean depth341 pdept(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1) ! bound by the ocean depth 282 342 ELSE 283 hd20(ji,jj) = 0._wp343 pdept(ji,jj) = 0._wp 284 344 ENDIF 285 !286 iid = ik28(ji,jj)287 IF( iid /= 1 ) THEN288 zztmp = gdept(ji,jj,iid ,Kmm) & ! linear interpolation289 & + ( gdept(ji,jj,iid+1,Kmm) - gdept(ji,jj,iid,Kmm) ) &290 & * ( 28.*tmask(ji,jj,iid+1) - ts(ji,jj,iid,jp_tem,Kmm) ) &291 & / ( ts(ji,jj,iid+1,jp_tem,Kmm) - ts(ji,jj,iid,jp_tem,Kmm) + (1.-tmask(ji,jj,1)) )292 hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1) ! bound by the ocean depth293 ELSE294 hd28(ji,jj) = 0._wp295 ENDIF296 297 345 END DO 298 346 END DO 299 CALL iom_put( "20d", hd20 ) ! depth of the 20 isotherm 300 CALL iom_put( "28d", hd28 ) ! depth of the 28 isotherm 301 302 ! ----------------------------- ! 303 ! Heat content of first 300 m ! 304 ! ----------------------------- ! 305 306 ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...) 307 ilevel = 0 308 zthick_0 = 0._wp 309 DO jk = 1, jpkm1 310 zthick_0 = zthick_0 + e3t_1d(jk) 311 IF( zthick_0 < 300. ) ilevel = jk 312 END DO 347 ! 348 END SUBROUTINE dia_hth_dep 349 350 351 SUBROUTINE dia_hth_htc( Kmm, pdep, pt, phtc ) 352 ! 353 INTEGER , INTENT(in) :: Kmm ! ocean time level index 354 REAL(wp), INTENT(in) :: pdep ! depth over the heat content 355 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pt 356 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phtc 357 ! 358 INTEGER :: ji, jj, jk, ik 359 REAL(wp), DIMENSION(jpi,jpj) :: zthick 360 INTEGER , DIMENSION(jpi,jpj) :: ilevel 361 362 313 363 ! surface boundary condition 314 IF( ln_linssh ) THEN ; zthick(:,:) = ssh(:,:,Kmm) ; htc3(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) * tmask(:,:,1) 315 ELSE ; zthick(:,:) = 0._wp ; htc3(:,:) = 0._wp 364 365 IF( .NOT. ln_linssh ) THEN ; zthick(:,:) = 0._wp ; phtc(:,:) = 0._wp 366 ELSE ; zthick(:,:) = ssh(:,:,Kmm) ; phtc(:,:) = pt(:,:,1) * ssh(:,:,Kmm) * tmask(:,:,1) 316 367 ENDIF 317 ! integration down to ilevel 318 DO jk = 1, ilevel 319 zthick(:,:) = zthick(:,:) + e3t(:,:,jk,Kmm) 320 htc3 (:,:) = htc3 (:,:) + e3t(:,:,jk,Kmm) * ts(:,:,jk,jp_tem,Kmm) * tmask(:,:,jk) 321 END DO 322 ! deepest layer 323 zthick(:,:) = 300. - zthick(:,:) ! remaining thickness to reach 300m 368 ! 369 ilevel(:,:) = 1 370 DO jk = 2, jpkm1 371 DO jj = 1, jpj 372 DO ji = 1, jpi 373 IF( ( gdept(ji,jj,jk,Kmm) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN 374 ilevel(ji,jj) = jk 375 zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm) 376 phtc (ji,jj) = phtc (ji,jj) + e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk) 377 ENDIF 378 ENDDO 379 ENDDO 380 ENDDO 381 ! 324 382 DO jj = 1, jpj 325 383 DO ji = 1, jpi 326 htc3(ji,jj) = htc3(ji,jj) + ts(ji,jj,ilevel+1,jp_tem,Kmm) & 327 & * MIN( e3t(ji,jj,ilevel+1,Kmm), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1) 384 ik = ilevel(ji,jj) 385 zthick(ji,jj) = pdep - zthick(ji,jj) ! remaining thickness to reach depht pdep 386 phtc(ji,jj) = phtc(ji,jj) + pt(ji,jj,ik+1) * MIN( e3t(ji,jj,ik+1,Kmm), zthick(ji,jj) ) & 387 * tmask(ji,jj,ik+1) 328 388 END DO 329 END DO 330 ! from temperature to heat contain 331 zcoef = rau0 * rcp 332 htc3(:,:) = zcoef * htc3(:,:) 333 CALL iom_put( "hc300", htc3 ) ! first 300m heat content 334 ! 335 IF( ln_timing ) CALL timing_stop('dia_hth') 336 ! 337 END SUBROUTINE dia_hth 338 339 #else 340 !!---------------------------------------------------------------------- 341 !! Default option : Empty module 342 !!---------------------------------------------------------------------- 343 LOGICAL , PUBLIC, PARAMETER :: lk_diahth = .FALSE. !: thermocline-20d depths flag 344 CONTAINS 345 SUBROUTINE dia_hth( kt, Kmm ) ! Empty routine 346 IMPLICIT NONE 347 INTEGER, INTENT( in ) :: kt 348 INTEGER, INTENT( in ) :: Kmm 349 WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt 350 END SUBROUTINE dia_hth 351 #endif 389 ENDDO 390 ! 391 ! 392 END SUBROUTINE dia_hth_htc 352 393 353 394 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
