Changeset 1577 for trunk/NEMO/OPA_SRC/DIA/diahth.F90
- Timestamp:
- 2009-08-04T16:56:59+02:00 (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DIA/diahth.F90
r1551 r1577 9 9 !! ! 1999-07 (E. Guilyardi) hd28 + heat content 10 10 !! 8.5 ! 2002-06 (G. Madec) F90: Free form and module 11 !! NEMO 3.2 ! 2009-07 (S. Masson) hc300 bugfix + cleaning 11 !! NEMO 3.2 ! 2009-07 (S. Masson) hc300 bugfix + cleaning + add new diag 12 12 !!---------------------------------------------------------------------- 13 13 … … 16 16 !! 'key_diahth' : thermocline depth diag. 17 17 !!---------------------------------------------------------------------- 18 !! dia_hth : Compute diagnostics associated with the thermocline18 !! dia_hth : Compute varius diagnostics associated with the mixed layer 19 19 !!---------------------------------------------------------------------- 20 20 !! * Modules used 21 21 USE oce ! ocean dynamics and tracers 22 22 USE dom_oce ! ocean space and time domain 23 USE zdf_oce ! ocean vertical physics 23 24 USE phycst ! physical constants 24 25 USE in_out_manager ! I/O manager … … 32 33 33 34 !! * Shared module variables 34 LOGICAL , PUBLIC, PARAMETER :: lk_diahth = .TRUE. !: thermocline-20d depths flag 35 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 36 hth , & !: depth of the max vertical temperature gradient (m) 37 hd20 , & !: depth of 20 C isotherm (m) 38 hd28 , & !: depth of 28 C isotherm (m) 39 htc3 !: heat content of first 300 m 35 LOGICAL , PUBLIC, PARAMETER :: lk_diahth = .TRUE. !: thermocline-20d depths flag 36 ! note: following variables should move to local variables once iom_put is always used 37 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hmld !: mixing layer depth (turbocline) [m] 38 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hth !: depth of the max vertical temperature gradient [m] 39 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hd20 !: depth of 20 C isotherm [m] 40 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hd28 !: depth of 28 C isotherm [m] 41 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: htc3 !: heat content of first 300 m [W] 40 42 41 43 !! * Substitutions … … 53 55 !! *** ROUTINE dia_hth *** 54 56 !! 55 !! ** Purpose : 56 !! Computes the depth of strongest vertical temperature gradient 57 !! Computes the depth of the 20 degree isotherm 58 !! Computes the depth of the 28 degree isotherm 59 !! Computes the heat content of first 300 m 57 !! ** Purpose : Computes 58 !! the mixing layer depth (turbocline): avt = 5.e-4 59 !! the depth of strongest vertical temperature gradient 60 !! the mixed layer depth with density criteria: rho = rho(10m or surf) + 0.03(or 0.01) 61 !! the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2 62 !! the top of the thermochine: tn = tn(10m) - ztem2 63 !! the pycnocline depth with density criteria equivalent to a temperature variation 64 !! rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 65 !! the barrier layer thickness 66 !! the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) ) 67 !! the depth of the 20 degree isotherm (linear interpolation) 68 !! the depth of the 28 degree isotherm (linear interpolation) 69 !! the heat content of first 300 m 60 70 !! 61 71 !! ** Method : … … 66 76 INTEGER :: ji, jj, jk ! dummy loop arguments 67 77 INTEGER :: iid, iif, ilevel ! temporary integers 68 INTEGER, DIMENSION(jpi,jpj) :: ikc ! levels 69 REAL(wp) :: zd, zthick_0, zcoef ! temporary scalars 70 REAL(wp), DIMENSION(jpi,jpj) :: zthick 71 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdzt 78 INTEGER, DIMENSION(jpi,jpj) :: ik20, ik28 ! levels 79 REAL(wp) :: zavt5 = 5.e-4_wp ! Kz criterion for the turbocline depth 80 REAL(wp) :: zrho3 = 0.03_wp ! density criterion for mixed layer depth 81 REAL(wp) :: zrho1 = 0.01_wp ! density criterion for mixed layer depth 82 REAL(wp) :: ztem2 = 0.2_wp ! temperature criterion for mixed layer depth 83 REAL(wp) :: zthick_0, zcoef ! temporary scalars 84 REAL(wp) :: zztmp, zzdep ! temporary scalars inside do loop 85 REAL(wp) :: zu, zv, zw, zut, zvt ! temporary workspace 86 REAL(wp), DIMENSION(jpi,jpj) :: zabs2 ! MLD: abs( tn - tn(10m) ) = ztem2 87 REAL(wp), DIMENSION(jpi,jpj) :: ztm2 ! Top of thermocline: tn = tn(10m) - ztem2 88 REAL(wp), DIMENSION(jpi,jpj) :: zrho10_3 ! MLD: rho = rho10m + zrho3 89 REAL(wp), DIMENSION(jpi,jpj) :: zpycn ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 90 REAL(wp), DIMENSION(jpi,jpj) :: ztinv ! max of temperature inversion 91 REAL(wp), DIMENSION(jpi,jpj) :: zdepinv ! depth of temperature inversion 92 REAL(wp), DIMENSION(jpi,jpj) :: zrho0_3 ! MLD rho = rho(surf) = 0.03 93 REAL(wp), DIMENSION(jpi,jpj) :: zrho0_1 ! MLD rho = rho(surf) = 0.01 94 REAL(wp), DIMENSION(jpi,jpj) :: zmaxdzT ! max of dT/dz 95 REAL(wp), DIMENSION(jpi,jpj) :: zthick ! vertical integration thickness 96 REAL(wp), DIMENSION(jpi,jpj) :: zdelr ! delta rho equivalent to deltaT = 0.2 72 97 !!---------------------------------------------------------------------- 73 98 … … 79 104 ENDIF 80 105 81 ! -------------------------- ! 82 ! Depth of the thermocline ! 83 ! -------------------------- ! 84 ! The depth of the thermocline is defined as the depth of the 85 ! strongest vertical temperature gradient 86 zdzt(:,:,1) = 0.e0 87 DO jk = 2, jpk ! vertical gradient of temperature 88 zdzt(:,:,jk) = ( tn(:,:,jk-1) - tn(:,:,jk) ) / fse3w(:,:,jk) * tmask(:,:,jk) 89 END DO 106 ! initialization 107 ztinv (:,:) = 0.e0_wp 108 zdepinv(:,:) = 0.e0_wp 109 zmaxdzT(:,:) = 0.e0_wp 90 110 DO jj = 1, jpj 91 111 DO ji = 1, jpi 92 ilevel = MAXLOC( zdzt(ji,jj,:), dim= 1 ) ! level of maximum vertical temperature gradient 93 hth(ji,jj) = fsdepw(ji,jj,ilevel) ! depth of the thermocline 94 END DO 95 END DO 96 97 CALL iom_put( "thermod", hth ) ! depth of the thermocline 98 99 ! ----------------------- ! 100 ! Depth of 20C isotherm ! 101 ! ----------------------- ! 112 zztmp = bathy(ji,jj) 113 hth (ji,jj) = zztmp 114 hmld (ji,jj) = zztmp 115 zabs2 (ji,jj) = zztmp 116 ztm2 (ji,jj) = zztmp 117 zrho10_3(ji,jj) = zztmp 118 zpycn (ji,jj) = zztmp 119 END DO 120 END DO 121 IF( nla10 > 1 ) THEN 122 DO jj = 1, jpj 123 DO ji = 1, jpi 124 zztmp = bathy(ji,jj) 125 zrho0_3(ji,jj) = zztmp 126 zrho0_1(ji,jj) = zztmp 127 END DO 128 END DO 129 ENDIF 102 130 103 ! search last level above 20C (beware temperature is not always decreasing with depth) 104 ikc(:,:) = MAXLOC( fsdept(:,:,1:jpkm1), mask = tn(:,:,1:jpkm1) >= 20., dim = 3 ) ! ikc = 0 if mask always equal to false 105 ! Depth of 20C isotherm, linear interpolation 131 ! Preliminary computation 132 ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 133 DO jj=1, jpj 134 DO ji=1, jpi 135 IF( tmask(ji,jj,nla10) == 1. ) THEN 136 zu = 1779.50 + 11.250*tn(ji,jj,nla10) - 3.80*sn(ji,jj,nla10) - 0.0745*tn(ji,jj,nla10)*tn(ji,jj,nla10) & 137 & - 0.0100*tn(ji,jj,nla10)*sn(ji,jj,nla10) 138 zv = 5891.00 + 38.000*tn(ji,jj,nla10) + 3.00*sn(ji,jj,nla10) - 0.3750*tn(ji,jj,nla10)*tn(ji,jj,nla10) 139 zut = 11.25 - 0.149*tn(ji,jj,nla10) - 0.01*sn(ji,jj,nla10) 140 zvt = 38.00 - 0.750*tn(ji,jj,nla10) 141 zw = (zu + 0.698*zv) * (zu + 0.698*zv) 142 zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 143 ELSE 144 zdelr(ji,jj) = 0.e0 145 ENDIF 146 END DO 147 END DO 148 149 ! ------------------------------------------------------------- ! 150 ! thermocline depth: strongest vertical gradient of temperature ! 151 ! turbocline depth (mixing layer depth): avt = zavt5 ! 152 ! MLD: rho = rho(1) + zrho3 ! 153 ! MLD: rho = rho(1) + zrho1 ! 154 ! ------------------------------------------------------------- ! 155 DO jk = jpkm1, 2, -1 ! loop from bottom to 2 156 DO jj = 1, jpj 157 DO ji = 1, jpi 158 159 zzdep = fsdepw(ji,jj,jk) 160 zztmp = ( tn(ji,jj,jk-1) - tn(ji,jj,jk) ) / zzdep * tmask(ji,jj,jk) ! vertical gradient of temperature (dT/dz) 161 zzdep = zzdep * tmask(ji,jj,1) 162 163 IF( zztmp > zmaxdzT(ji,jj) ) THEN 164 zmaxdzT(ji,jj) = zztmp ; hth (ji,jj) = zzdep ! max and depth of dT/dz 165 ENDIF 166 167 IF( avt (ji,jj,jk) < zavt5 ) hmld (ji,jj) = zzdep ! avt < zavt5 168 169 IF( nla10 > 1 ) THEN 170 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1) ! delta rho(1) 171 IF( zztmp > zrho3 ) zrho0_3(ji,jj) = zzdep ! > 0.03 172 IF( zztmp > zrho1 ) zrho0_1(ji,jj) = zzdep ! > 0.01 173 ENDIF 174 175 END DO 176 END DO 177 END DO 178 179 CALL iom_put( "mlddzt", hth ) ! depth of the thermocline 180 CALL iom_put( "mldkz5", hmld ) ! turbocline depth 181 IF( nla10 > 1 ) THEN 182 CALL iom_put( "mldr0_3", zrho0_3 ) ! MLD delta rho(surf) = 0.03 183 CALL iom_put( "mldr0_1", zrho0_1 ) ! MLD delta rho(surf) = 0.01 184 ENDIF 185 186 ! ------------------------------------------------------------- ! 187 ! MLD: abs( tn - tn(10m) ) = ztem2 ! 188 ! Top of thermocline: tn = tn(10m) - ztem2 ! 189 ! MLD: rho = rho10m + zrho3 ! 190 ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) ! 191 ! temperature inversion: max( 0, max of tn - tn(10m) ) ! 192 ! depth of temperature inversion ! 193 ! ------------------------------------------------------------- ! 194 DO jk = jpkm1, nlb10, -1 ! loop from bottom to nlb10 195 DO jj = 1, jpj 196 DO ji = 1, jpi 197 198 zzdep = fsdepw(ji,jj,jk) * tmask(ji,jj,1) 199 200 zztmp = tn(ji,jj,nla10) - tn(ji,jj,jk) ! - delta T(10m) 201 IF( ABS(zztmp) > ztem2 ) zabs2 (ji,jj) = zzdep ! abs > 0.2 202 IF( zztmp > ztem2 ) ztm2 (ji,jj) = zzdep ! > 0.2 203 zztmp = -zztmp ! delta T(10m) 204 IF( zztmp > ztinv(ji,jj) ) THEN ! temperature inversion 205 ztinv(ji,jj) = zztmp ; zdepinv (ji,jj) = zzdep ! max value and depth 206 ENDIF 207 208 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10) ! delta rho(10m) 209 IF( zztmp > zrho3 ) zrho10_3(ji,jj) = zzdep ! > 0.03 210 IF( zztmp > zdelr(ji,jj) ) zpycn (ji,jj) = zzdep ! > equi. delta T(10m) - 0.2 211 212 END DO 213 END DO 214 END DO 215 216 CALL iom_put( "mld|dt|" , zabs2 ) ! MLD abs(delta t) - 0.2 217 CALL iom_put( "topthdep", ztm2 ) ! T(10) - 0.2 218 CALL iom_put( "mldr10_3", zrho10_3 ) ! MLD delta rho(10m) = 0.03 219 CALL iom_put( "pycndep" , zpycn ) ! MLD delta rho equi. delta T(10m) = 0.2 220 CALL iom_put( "BLT" , ztm2 - zpycn ) ! Barrier Layer Thickness 221 CALL iom_put( "tinv" , ztinv ) ! max. temp. inv. (t10 ref) 222 CALL iom_put( "depti" , zdepinv ) ! depth of max. temp. inv. (t10 ref) 223 224 225 ! ----------------------------------- ! 226 ! search deepest level above 20C/28C ! 227 ! ----------------------------------- ! 228 ik20(:,:) = 1 229 ik28(:,:) = 1 230 DO jk = 1, jpkm1 ! beware temperature is not always decreasing with depth => loop from top to bottom 231 DO jj = 1, jpj 232 DO ji = 1, jpi 233 zztmp = tn(ji,jj,jk) 234 IF( zztmp >= 20. ) ik20(ji,jj) = jk 235 IF( zztmp >= 28. ) ik28(ji,jj) = jk 236 END DO 237 END DO 238 END DO 239 240 ! --------------------------- ! 241 ! Depth of 20C/28C isotherm ! 242 ! --------------------------- ! 106 243 DO jj = 1, jpj 107 244 DO ji = 1, jpi 108 iid = MAX( 1, ikc(ji,jj) ) ! make sure that iid /= 0 (these points are masked later) 109 zd = fsdept(ji,jj,iid) + ( fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid) ) & 110 & * ( 20.*tmask(ji,jj,iid+1) - tn(ji,jj,iid) ) & 111 & / ( tn(ji,jj,iid+1) - tn(ji,jj,iid) + ( 1. - tmask(ji,jj,1) ) ) 112 ! bound by the ocean depth, minimum value, first T-point depth 245 113 246 iif = mbathy(ji,jj) 114 hd20(ji,jj) = MIN( zd*tmask(ji,jj,1), fsdepw(ji,jj,iif) ) 115 END DO 116 END DO 117 WHERE(ikc == 0 ) hd20 = 0.e0 ! set to 0 points where tn is always < 20 (not very good for temporal mean...) 118 CALL iom_put( "20d", hd20 ) ! depth of the 20 isotherm 119 120 ! ----------------------- ! 121 ! Depth of 28C isotherm ! 122 ! ----------------------- ! 123 124 ! search last level above 28C (beware temperature is not always decreasing with depth) 125 ikc(:,:) = MAXLOC( fsdept(:,:,1:jpkm1), mask = tn(:,:,1:jpkm1) >= 28., dim = 3 ) ! ikc = 0 if mask always equal to false 126 ! Depth of 28C isotherm, linear interpolation 127 DO jj = 1, jpj 128 DO ji = 1, jpi 129 iid = MAX( 1, ikc(ji,jj) ) ! make sure that iid /= 0 (these points are masked later) 130 zd = fsdept(ji,jj,iid) + ( fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid) ) & 131 & * ( 28.*tmask(ji,jj,iid+1) - tn(ji,jj,iid) ) & 132 & / ( tn(ji,jj,iid+1) - tn(ji,jj,iid) + ( 1. - tmask(ji,jj,1) ) ) 133 ! bound by the ocean depth, minimum value, first T-point depth 134 iif = mbathy(ji,jj) 135 hd28(ji,jj) = MIN( zd*tmask(ji,jj,1), fsdepw(ji,jj,iif) ) 136 END DO 137 END DO 138 WHERE(ikc == 0 ) hd28 = 0.e0 ! set to 0 points where tn is always < 28 (not very good for temporal mean...) 139 CALL iom_put( "28d", hd28 ) ! depth of the 28 isotherm 247 zzdep = fsdepw(ji,jj,iif) 248 249 iid = ik20(ji,jj) 250 IF( iid /= 1 ) THEN 251 ! linear interpolation 252 zztmp = fsdept(ji,jj,iid ) & 253 & + ( fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid) ) & 254 & * ( 20.*tmask(ji,jj,iid+1) - tn(ji,jj,iid) ) & 255 & / ( tn(ji,jj,iid+1) - tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) ) 256 ! bound by the ocean depth, minimum value, first T-point depth 257 hd20(ji,jj) = MIN( zztmp*tmask(ji,jj,1), zzdep) 258 ELSE 259 hd20(ji,jj)=0. 260 ENDIF 261 262 iid = ik28(ji,jj) 263 IF( iid /= 1 ) THEN 264 ! linear interpolation 265 zztmp = fsdept(ji,jj,iid ) & 266 & + ( fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid) ) & 267 & * ( 28.*tmask(ji,jj,iid+1) - tn(ji,jj,iid) ) & 268 & / ( tn(ji,jj,iid+1) - tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) ) 269 ! bound by the ocean depth, minimum value, first T-point depth 270 hd28(ji,jj) = MIN( zztmp*tmask(ji,jj,1), zzdep ) 271 ELSE 272 hd28(ji,jj) = 0. 273 ENDIF 274 275 END DO 276 END DO 277 CALL iom_put( "20d", hd20 ) ! depth of the 20 isotherm 278 CALL iom_put( "28d", hd28 ) ! depth of the 28 isotherm 140 279 141 280 ! ----------------------------- ! … … 145 284 ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_0 to do this search...) 146 285 ilevel = 0 147 zthick_0 = 0.e0 286 zthick_0 = 0.e0_wp 148 287 DO jk = 1, jpkm1 149 288 zthick_0 = zthick_0 + e3t_0(jk) … … 151 290 END DO 152 291 ! surface boundary condition 153 IF( lk_vvl ) THEN ; zthick(:,:) = 0.e0 ; htc3(:,:) = 0.e0292 IF( lk_vvl ) THEN ; zthick(:,:) = 0.e0_wp ; htc3(:,:) = 0.e0_wp 154 293 ELSE ; zthick(:,:) = sshn(:,:) ; htc3(:,:) = tn(:,:,jk) * sshn(:,:) * tmask(:,:,jk) 155 294 ENDIF … … 169 308 zcoef = rau0 * rcp 170 309 htc3(:,:) = zcoef * htc3(:,:) 171 CALL iom_put( "hc300", htc3 ) ! first 300m hea at content310 CALL iom_put( "hc300", htc3 ) ! first 300m heat content 172 311 173 312
Note: See TracChangeset
for help on using the changeset viewer.