Changeset 12517
- Timestamp:
- 2020-03-05T16:11:56+01:00 (3 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE
- Files:
-
- 3 edited
-
C1D/step_c1d.F90 (modified) (2 diffs)
-
DIA/diahth.F90 (modified) (4 diffs)
-
step.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/C1D/step_c1d.F90
r11715 r12517 58 58 59 59 indic = 0 ! reset to no error condition 60 IF( kstp == nit000 ) CALL iom_init( "nemo") ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS) 60 IF( kstp == nit000 ) THEN 61 CALL iom_init( "nemo") ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS) 62 CALL dia_hth_init ! extra ML depth diagnostics, thermocline depths 63 END IF 61 64 IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init) 62 65 CALL iom_setkt( kstp - nit000 + 1, "nemo" ) ! say to iom that we are at time step kstp … … 86 89 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 87 90 CALL dia_wri( kstp ) ! ocean model: outputs 88 IF( l k_diahth ) CALL dia_hth( kstp ) ! Thermocline depth (20°C)91 IF( ll_diahth ) CALL dia_hth( kstp ) ! Thermocline depth (20°C) 89 92 90 93 -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/DIA/diahth.F90
r11715 r12517 5 5 !!====================================================================== 6 6 !! History : OPA ! 1994-09 (J.-P. Boulanger) Original code 7 !! ! 1996-11 (E. Guilyardi) OPA8 7 !! ! 1996-11 (E. Guilyardi) OPA8 8 8 !! ! 1997-08 (G. Madec) optimization 9 !! ! 1999-07 (E. Guilyardi) hd28 + heat content 9 !! ! 1999-07 (E. Guilyardi) hd28 + heat content 10 10 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 11 11 !! 3.2 ! 2009-07 (S. Masson) hc300 bugfix + cleaning + add new diag 12 !!---------------------------------------------------------------------- 13 #if defined key_diahth 14 !!---------------------------------------------------------------------- 15 !! 'key_diahth' : thermocline depth diag. 12 !! 4.0.1! 2020-01 (G. Nurser) remove need for key lk_diahth 16 13 !!---------------------------------------------------------------------- 17 14 !! dia_hth : Compute varius diagnostics associated with the mixed layer … … 24 21 USE lib_mpp ! MPP library 25 22 USE iom ! I/O library 26 USE timing ! preformance summary27 23 28 24 IMPLICIT NONE … … 30 26 31 27 PUBLIC dia_hth ! routine called by step.F90 32 PUBLIC dia_hth_alloc ! routine called by nemogcm.F90 33 34 LOGICAL , PUBLIC, PARAMETER :: lk_diahth = .TRUE. !: thermocline-20d depths flag 35 36 ! note: following variables should move to local variables once iom_put is always used 37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hth !: depth of the max vertical temperature gradient [m] 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hd20 !: depth of 20 C isotherm [m] 39 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hd28 !: depth of 28 C isotherm [m] 40 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: htc3 !: heat content of first 300 m [W] 28 PUBLIC dia_hth_init ! routine called by nemogcm.F90 29 30 LOGICAL, PUBLIC :: ll_diahth !: Compute further diagnostics of ML and thermocline depth 41 31 42 32 !!---------------------------------------------------------------------- … … 47 37 CONTAINS 48 38 49 FUNCTION dia_hth_alloc() 50 !!--------------------------------------------------------------------- 51 INTEGER :: dia_hth_alloc 52 !!--------------------------------------------------------------------- 39 40 SUBROUTINE dia_hth( kt ) 41 !!--------------------------------------------------------------------- 42 !! *** ROUTINE dia_hth *** 43 !! 44 !! ** Purpose : Computes 45 !! the mixing layer depth (turbocline): avt = 5.e-4 46 !! the depth of strongest vertical temperature gradient 47 !! the mixed layer depth with density criteria: rho = rho(10m or surf) + 0.03(or 0.01) 48 !! the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2 49 !! the top of the thermochine: tn = tn(10m) - ztem2 50 !! the pycnocline depth with density criteria equivalent to a temperature variation 51 !! rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 52 !! the barrier layer thickness 53 !! the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) ) 54 !! the depth of the 20 degree isotherm (linear interpolation) 55 !! the depth of the 28 degree isotherm (linear interpolation) 56 !! the heat content of first 300 m 57 !! 58 !! ** Method : 59 !!------------------------------------------------------------------- 60 INTEGER, INTENT( in ) :: kt ! ocean time-step index 61 !! 62 INTEGER :: ji, jj, jk ! dummy loop arguments 63 INTEGER :: iid, ilevel ! temporary integers 64 INTEGER, DIMENSION(jpi,jpj) :: ik20, ik28 ! levels 65 REAL(wp) :: zavt5 = 5.e-4_wp ! Kz criterion for the turbocline depth 66 REAL(wp) :: zrho3 = 0.03_wp ! density criterion for mixed layer depth 67 REAL(wp) :: zrho1 = 0.01_wp ! density criterion for mixed layer depth 68 REAL(wp) :: ztem2 = 0.2_wp ! temperature criterion for mixed layer depth 69 REAL(wp) :: zthick_0, zcoef ! temporary scalars 70 REAL(wp) :: zztmp, zzdep ! temporary scalars inside do loop 71 REAL(wp) :: zu, zv, zw, zut, zvt ! temporary workspace 72 REAL(wp), DIMENSION(jpi,jpj) :: zabs2 ! MLD: abs( tn - tn(10m) ) = ztem2 73 REAL(wp), DIMENSION(jpi,jpj) :: ztm2 ! Top of thermocline: tn = tn(10m) - ztem2 74 REAL(wp), DIMENSION(jpi,jpj) :: zrho10_3 ! MLD: rho = rho10m + zrho3 75 REAL(wp), DIMENSION(jpi,jpj) :: zpycn ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 76 REAL(wp), DIMENSION(jpi,jpj) :: ztinv ! max of temperature inversion 77 REAL(wp), DIMENSION(jpi,jpj) :: zdepinv ! depth of temperature inversion 78 REAL(wp), DIMENSION(jpi,jpj) :: zrho0_3 ! MLD rho = rho(surf) = 0.03 79 REAL(wp), DIMENSION(jpi,jpj) :: zrho0_1 ! MLD rho = rho(surf) = 0.01 80 REAL(wp), DIMENSION(jpi,jpj) :: zmaxdzT ! max of dT/dz 81 REAL(wp), DIMENSION(jpi,jpj) :: zthick ! vertical integration thickness 82 REAL(wp), DIMENSION(jpi,jpj) :: zdelr ! delta rho equivalent to deltaT = 0.2 83 ! note: following variables should move to local variables once iom_put is always used 84 REAL(wp), DIMENSION(jpi,jpj) :: zhth !: depth of the max vertical temperature gradient [m] 85 REAL(wp), DIMENSION(jpi,jpj) :: zhd20 !: depth of 20 C isotherm [m] 86 REAL(wp), DIMENSION(jpi,jpj) :: zhd28 !: depth of 28 C isotherm [m] 87 REAL(wp), DIMENSION(jpi,jpj) :: zhtc3 !: heat content of first 300 m [W] 88 89 IF (iom_use("mlddzt") .OR. iom_use("mldr0_3") .OR. iom_use("mldr0_1")) THEN 90 ! ------------------------------------------------------------- ! 91 ! thermocline depth: strongest vertical gradient of temperature ! 92 ! turbocline depth (mixing layer depth): avt = zavt5 ! 93 ! MLD: rho = rho(1) + zrho3 ! 94 ! MLD: rho = rho(1) + zrho1 ! 95 ! ------------------------------------------------------------- ! 96 zmaxdzT(:,:) = 0._wp 97 IF( nla10 > 1 ) THEN 98 DO jj = 1, jpj 99 DO ji = 1, jpi 100 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 101 zrho0_3(ji,jj) = zztmp 102 zrho0_1(ji,jj) = zztmp 103 zhth(ji,jj) = zztmp 104 END DO 105 END DO 106 ELSE IF (iom_use("mlddzt")) THEN 107 DO jj = 1, jpj 108 DO ji = 1, jpi 109 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 110 zhth(ji,jj) = zztmp 111 END DO 112 END DO 113 ELSE 114 zhth(:,:) = 0._wp 115 116 ENDIF 117 118 DO jk = jpkm1, 2, -1 ! loop from bottom to 2 119 DO jj = 1, jpj 120 DO ji = 1, jpi 121 ! 122 zzdep = gdepw_n(ji,jj,jk) 123 zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk) ! vertical gradient of temperature (dT/dz) 124 zzdep = zzdep * tmask(ji,jj,1) 125 126 IF( zztmp > zmaxdzT(ji,jj) ) THEN 127 zmaxdzT(ji,jj) = zztmp ; zhth (ji,jj) = zzdep ! max and depth of dT/dz 128 ENDIF 129 130 IF( nla10 > 1 ) THEN 131 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1) ! delta rho(1) 132 IF( zztmp > zrho3 ) zrho0_3(ji,jj) = zzdep ! > 0.03 133 IF( zztmp > zrho1 ) zrho0_1(ji,jj) = zzdep ! > 0.01 134 ENDIF 135 136 END DO 137 END DO 138 END DO 139 140 IF (iom_use("mlddzt")) CALL iom_put( "mlddzt", zhth*tmask(:,:,1) ) ! depth of the thermocline 141 IF( nla10 > 1 ) THEN 142 IF (iom_use("mldr0_3")) CALL iom_put( "mldr0_3", zrho0_3*tmask(:,:,1) ) ! MLD delta rho(surf) = 0.03 143 IF (iom_use("mldr0_1")) CALL iom_put( "mldr0_1", zrho0_1*tmask(:,:,1) ) ! MLD delta rho(surf) = 0.01 144 ENDIF 145 ENDIF 146 147 IF (iom_use("mld_dt02") .OR. iom_use("topthdep") .OR. iom_use("mldr10_3") .OR. & 148 & iom_use("pycndep") .OR. iom_use("tinv") .OR. iom_use("depti")) THEN 149 DO jj = 1, jpj 150 DO ji = 1, jpi 151 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 152 zabs2 (ji,jj) = zztmp 153 ztm2 (ji,jj) = zztmp 154 zrho10_3(ji,jj) = zztmp 155 zpycn (ji,jj) = zztmp 156 END DO 157 END DO 158 ztinv (:,:) = 0._wp 159 zdepinv(:,:) = 0._wp 160 161 IF (iom_use("pycndep")) THEN 162 ! Preliminary computation 163 ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 164 DO jj = 1, jpj 165 DO ji = 1, jpi 166 IF( tmask(ji,jj,nla10) == 1. ) THEN 167 zu = 1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80 * tsn(ji,jj,nla10,jp_sal) & 168 & - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) & 169 & - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal) 170 zv = 5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00 * tsn(ji,jj,nla10,jp_sal) & 171 & - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) 172 zut = 11.25 - 0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01 * tsn(ji,jj,nla10,jp_sal) 173 zvt = 38.00 - 0.750 * tsn(ji,jj,nla10,jp_tem) 174 zw = (zu + 0.698*zv) * (zu + 0.698*zv) 175 zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 176 ELSE 177 zdelr(ji,jj) = 0._wp 178 ENDIF 179 END DO 180 END DO 181 ELSE 182 zdelr(:,:) = 0._wp 183 ENDIF 184 185 ! ------------------------------------------------------------- ! 186 ! MLD: abs( tn - tn(10m) ) = ztem2 ! 187 ! Top of thermocline: tn = tn(10m) - ztem2 ! 188 ! MLD: rho = rho10m + zrho3 ! 189 ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) ! 190 ! temperature inversion: max( 0, max of tn - tn(10m) ) ! 191 ! depth of temperature inversion ! 192 ! ------------------------------------------------------------- ! 193 DO jk = jpkm1, nlb10, -1 ! loop from bottom to nlb10 194 DO jj = 1, jpj 195 DO ji = 1, jpi 196 ! 197 zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1) 198 ! 199 zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem) ! - delta T(10m) 200 IF( ABS(zztmp) > ztem2 ) zabs2 (ji,jj) = zzdep ! abs > 0.2 201 IF( zztmp > ztem2 ) ztm2 (ji,jj) = zzdep ! > 0.2 202 zztmp = -zztmp ! delta T(10m) 203 IF( zztmp > ztinv(ji,jj) ) THEN ! temperature inversion 204 ztinv(ji,jj) = zztmp ; zdepinv (ji,jj) = zzdep ! max value and depth 205 ENDIF 206 207 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10) ! delta rho(10m) 208 IF( zztmp > zrho3 ) zrho10_3(ji,jj) = zzdep ! > 0.03 209 IF( zztmp > zdelr(ji,jj) ) zpycn (ji,jj) = zzdep ! > equi. delta T(10m) - 0.2 210 ! 211 END DO 212 END DO 213 END DO 214 215 IF (iom_use("mld_dt02")) CALL iom_put( "mld_dt02", zabs2*tmask(:,:,1) ) ! MLD abs(delta t) - 0.2 216 IF (iom_use("topthdep")) CALL iom_put( "topthdep", ztm2*tmask(:,:,1) ) ! T(10) - 0.2 217 IF (iom_use("mldr10_3")) CALL iom_put( "mldr10_3", zrho10_3*tmask(:,:,1) ) ! MLD delta rho(10m) = 0.03 218 IF (iom_use("pycndep")) CALL iom_put( "pycndep" , zpycn*tmask(:,:,1) ) ! MLD delta rho equi. delta T(10m) = 0.2 219 IF (iom_use("tinv")) CALL iom_put( "tinv" , ztinv*tmask(:,:,1) ) ! max. temp. inv. (t10 ref) 220 IF (iom_use("depti")) CALL iom_put( "depti" , zdepinv*tmask(:,:,1) ) ! depth of max. temp. inv. (t10 ref) 221 ENDIF 222 223 IF(iom_use("20d") .OR. iom_use("28d")) THEN 224 ! ----------------------------------- ! 225 ! search deepest level above 20C/28C ! 226 ! ----------------------------------- ! 227 ik20(:,:) = 1 228 ik28(:,:) = 1 229 DO jk = 1, jpkm1 ! beware temperature is not always decreasing with depth => loop from top to bottom 230 DO jj = 1, jpj 231 DO ji = 1, jpi 232 zztmp = tsn(ji,jj,jk,jp_tem) 233 IF( zztmp >= 20. ) ik20(ji,jj) = jk 234 IF( zztmp >= 28. ) ik28(ji,jj) = jk 235 END DO 236 END DO 237 END DO 238 239 ! --------------------------- ! 240 ! Depth of 20C/28C isotherm ! 241 ! --------------------------- ! 242 DO jj = 1, jpj 243 DO ji = 1, jpi 244 ! 245 zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1) ! depth of the oean bottom 246 ! 247 iid = ik20(ji,jj) 248 IF( iid /= 1 ) THEN 249 zztmp = gdept_n(ji,jj,iid ) & ! linear interpolation 250 & + ( gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid) ) & 251 & * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem) ) & 252 & / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 253 zhd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1) ! bound by the ocean depth 254 ELSE 255 zhd20(ji,jj) = 0._wp 256 ENDIF 257 ! 258 iid = ik28(ji,jj) 259 IF( iid /= 1 ) THEN 260 zztmp = gdept_n(ji,jj,iid ) & ! linear interpolation 261 & + ( gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid) ) & 262 & * ( 28.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem) ) & 263 & / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 264 zhd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1) ! bound by the ocean depth 265 ELSE 266 zhd28(ji,jj) = 0._wp 267 ENDIF 268 269 END DO 270 END DO 271 CALL iom_put( "20d", zhd20 ) ! depth of the 20 isotherm 272 CALL iom_put( "28d", zhd28 ) ! depth of the 28 isotherm 273 ENDIF 274 275 ! ----------------------------- ! 276 ! Heat content of first 300 m ! 277 ! ----------------------------- ! 278 IF (iom_use("hc300")) THEN 279 ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...) 280 ilevel = 0 281 zthick_0 = 0._wp 282 DO jk = 1, jpkm1 283 zthick_0 = zthick_0 + e3t_1d(jk) 284 IF( zthick_0 < 300. ) ilevel = jk 285 END DO 286 ! surface boundary condition 287 IF( ln_linssh ) THEN ; zthick(:,:) = sshn(:,:) ; zhtc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 288 ELSE ; zthick(:,:) = 0._wp ; zhtc3(:,:) = 0._wp 289 ENDIF 290 ! integration down to ilevel 291 DO jk = 1, ilevel 292 zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk) 293 zhtc3 (:,:) = zhtc3 (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk) 294 END DO 295 ! deepest layer 296 zthick(:,:) = 300. - zthick(:,:) ! remaining thickness to reach 300m 297 DO jj = 1, jpj 298 DO ji = 1, jpi 299 zhtc3(ji,jj) = zhtc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) & 300 & * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1) 301 END DO 302 END DO 303 ! from temperature to heat contain 304 zcoef = rau0 * rcp 305 zhtc3(:,:) = zcoef * zhtc3(:,:) 306 CALL iom_put( "hc300", zhtc3*tmask(:,:,1) ) ! first 300m heat content 307 ENDIF 308 ! 309 END SUBROUTINE dia_hth 310 311 312 SUBROUTINE dia_hth_init 313 !!--------------------------------------------------------------------------- 314 !! *** ROUTINE dia_hth_init *** 315 !! 316 !! ** Purpose: Initialization for ML and thermocline depths 317 !! 318 !! ** Action : If any upper ocean diagnostic required by xml file, set in dia_hth 319 !!--------------------------------------------------------------------------- 320 ! 321 IF(lwp) THEN 322 WRITE(numout,*) 323 WRITE(numout,*) 'dia_hth_init : heat and salt budgets diagnostics' 324 WRITE(numout,*) '~~~~~~~~~~~~ ' 325 ENDIF 326 ll_diahth = iom_use("mlddzt") .OR. iom_use("mldr0_3") .OR. iom_use("mldr0_1") .OR. & 327 & iom_use("mld_dt02") .OR. iom_use("topthdep") .OR. iom_use("mldr10_3") .OR. & 328 & iom_use("pycndep") .OR. iom_use("tinv") .OR. iom_use("depti").OR. & 329 & iom_use("20d") .OR. iom_use("28d") .OR. iom_use("hc300") 330 IF(lwp) THEN 331 WRITE(numout,*) ' output upper ocean diagnostics (T) or not (F) ll_diahth = ', ll_diahth 332 ENDIF 53 333 ! 54 ALLOCATE( hth(jpi,jpj), hd20(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), STAT=dia_hth_alloc ) 55 ! 56 CALL mpp_sum ( 'diahth', dia_hth_alloc ) 57 IF(dia_hth_alloc /= 0) CALL ctl_stop( 'STOP', 'dia_hth_alloc: failed to allocate arrays.' ) 58 ! 59 END FUNCTION dia_hth_alloc 60 61 62 SUBROUTINE dia_hth( kt ) 63 !!--------------------------------------------------------------------- 64 !! *** ROUTINE dia_hth *** 65 !! 66 !! ** Purpose : Computes 67 !! the mixing layer depth (turbocline): avt = 5.e-4 68 !! the depth of strongest vertical temperature gradient 69 !! the mixed layer depth with density criteria: rho = rho(10m or surf) + 0.03(or 0.01) 70 !! the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2 71 !! the top of the thermochine: tn = tn(10m) - ztem2 72 !! the pycnocline depth with density criteria equivalent to a temperature variation 73 !! rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 74 !! the barrier layer thickness 75 !! the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) ) 76 !! the depth of the 20 degree isotherm (linear interpolation) 77 !! the depth of the 28 degree isotherm (linear interpolation) 78 !! the heat content of first 300 m 79 !! 80 !! ** Method : 81 !!------------------------------------------------------------------- 82 INTEGER, INTENT( in ) :: kt ! ocean time-step index 83 !! 84 INTEGER :: ji, jj, jk ! dummy loop arguments 85 INTEGER :: iid, ilevel ! temporary integers 86 INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ik20, ik28 ! levels 87 REAL(wp) :: zavt5 = 5.e-4_wp ! Kz criterion for the turbocline depth 88 REAL(wp) :: zrho3 = 0.03_wp ! density criterion for mixed layer depth 89 REAL(wp) :: zrho1 = 0.01_wp ! density criterion for mixed layer depth 90 REAL(wp) :: ztem2 = 0.2_wp ! temperature criterion for mixed layer depth 91 REAL(wp) :: zthick_0, zcoef ! temporary scalars 92 REAL(wp) :: zztmp, zzdep ! temporary scalars inside do loop 93 REAL(wp) :: zu, zv, zw, zut, zvt ! temporary workspace 94 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zabs2 ! MLD: abs( tn - tn(10m) ) = ztem2 95 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ztm2 ! Top of thermocline: tn = tn(10m) - ztem2 96 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zrho10_3 ! MLD: rho = rho10m + zrho3 97 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zpycn ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 98 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ztinv ! max of temperature inversion 99 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdepinv ! depth of temperature inversion 100 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zrho0_3 ! MLD rho = rho(surf) = 0.03 101 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zrho0_1 ! MLD rho = rho(surf) = 0.01 102 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zmaxdzT ! max of dT/dz 103 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zthick ! vertical integration thickness 104 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdelr ! delta rho equivalent to deltaT = 0.2 105 !!---------------------------------------------------------------------- 106 IF( ln_timing ) CALL timing_start('dia_hth') 107 108 IF( kt == nit000 ) THEN 109 ! ! allocate dia_hth array 110 IF( dia_hth_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard arrays' ) 111 112 IF(.NOT. ALLOCATED(ik20) ) THEN 113 ALLOCATE(ik20(jpi,jpj), ik28(jpi,jpj), & 114 & zabs2(jpi,jpj), & 115 & ztm2(jpi,jpj), & 116 & zrho10_3(jpi,jpj),& 117 & zpycn(jpi,jpj), & 118 & ztinv(jpi,jpj), & 119 & zdepinv(jpi,jpj), & 120 & zrho0_3(jpi,jpj), & 121 & zrho0_1(jpi,jpj), & 122 & zmaxdzT(jpi,jpj), & 123 & zthick(jpi,jpj), & 124 & zdelr(jpi,jpj), STAT=ji) 125 CALL mpp_sum('diahth', ji) 126 IF( ji /= 0 ) CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard ocean arrays' ) 127 END IF 128 129 IF(lwp) WRITE(numout,*) 130 IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth' 131 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 132 IF(lwp) WRITE(numout,*) 133 ENDIF 134 135 ! initialization 136 ztinv (:,:) = 0._wp 137 zdepinv(:,:) = 0._wp 138 zmaxdzT(:,:) = 0._wp 139 DO jj = 1, jpj 140 DO ji = 1, jpi 141 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 142 hth (ji,jj) = zztmp 143 zabs2 (ji,jj) = zztmp 144 ztm2 (ji,jj) = zztmp 145 zrho10_3(ji,jj) = zztmp 146 zpycn (ji,jj) = zztmp 147 END DO 148 END DO 149 IF( nla10 > 1 ) THEN 150 DO jj = 1, jpj 151 DO ji = 1, jpi 152 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 153 zrho0_3(ji,jj) = zztmp 154 zrho0_1(ji,jj) = zztmp 155 END DO 156 END DO 157 ENDIF 158 159 ! Preliminary computation 160 ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 161 DO jj = 1, jpj 162 DO ji = 1, jpi 163 IF( tmask(ji,jj,nla10) == 1. ) THEN 164 zu = 1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80 * tsn(ji,jj,nla10,jp_sal) & 165 & - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) & 166 & - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal) 167 zv = 5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00 * tsn(ji,jj,nla10,jp_sal) & 168 & - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) 169 zut = 11.25 - 0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01 * tsn(ji,jj,nla10,jp_sal) 170 zvt = 38.00 - 0.750 * tsn(ji,jj,nla10,jp_tem) 171 zw = (zu + 0.698*zv) * (zu + 0.698*zv) 172 zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 173 ELSE 174 zdelr(ji,jj) = 0._wp 175 ENDIF 176 END DO 177 END DO 178 179 ! ------------------------------------------------------------- ! 180 ! thermocline depth: strongest vertical gradient of temperature ! 181 ! turbocline depth (mixing layer depth): avt = zavt5 ! 182 ! MLD: rho = rho(1) + zrho3 ! 183 ! MLD: rho = rho(1) + zrho1 ! 184 ! ------------------------------------------------------------- ! 185 DO jk = jpkm1, 2, -1 ! loop from bottom to 2 186 DO jj = 1, jpj 187 DO ji = 1, jpi 188 ! 189 zzdep = gdepw_n(ji,jj,jk) 190 zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk) ! vertical gradient of temperature (dT/dz) 191 zzdep = zzdep * tmask(ji,jj,1) 192 193 IF( zztmp > zmaxdzT(ji,jj) ) THEN 194 zmaxdzT(ji,jj) = zztmp ; hth (ji,jj) = zzdep ! max and depth of dT/dz 195 ENDIF 196 197 IF( nla10 > 1 ) THEN 198 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1) ! delta rho(1) 199 IF( zztmp > zrho3 ) zrho0_3(ji,jj) = zzdep ! > 0.03 200 IF( zztmp > zrho1 ) zrho0_1(ji,jj) = zzdep ! > 0.01 201 ENDIF 202 203 END DO 204 END DO 205 END DO 206 207 CALL iom_put( "mlddzt", hth ) ! depth of the thermocline 208 IF( nla10 > 1 ) THEN 209 CALL iom_put( "mldr0_3", zrho0_3 ) ! MLD delta rho(surf) = 0.03 210 CALL iom_put( "mldr0_1", zrho0_1 ) ! MLD delta rho(surf) = 0.01 211 ENDIF 212 213 ! ------------------------------------------------------------- ! 214 ! MLD: abs( tn - tn(10m) ) = ztem2 ! 215 ! Top of thermocline: tn = tn(10m) - ztem2 ! 216 ! MLD: rho = rho10m + zrho3 ! 217 ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) ! 218 ! temperature inversion: max( 0, max of tn - tn(10m) ) ! 219 ! depth of temperature inversion ! 220 ! ------------------------------------------------------------- ! 221 DO jk = jpkm1, nlb10, -1 ! loop from bottom to nlb10 222 DO jj = 1, jpj 223 DO ji = 1, jpi 224 ! 225 zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1) 226 ! 227 zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem) ! - delta T(10m) 228 IF( ABS(zztmp) > ztem2 ) zabs2 (ji,jj) = zzdep ! abs > 0.2 229 IF( zztmp > ztem2 ) ztm2 (ji,jj) = zzdep ! > 0.2 230 zztmp = -zztmp ! delta T(10m) 231 IF( zztmp > ztinv(ji,jj) ) THEN ! temperature inversion 232 ztinv(ji,jj) = zztmp ; zdepinv (ji,jj) = zzdep ! max value and depth 233 ENDIF 234 235 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10) ! delta rho(10m) 236 IF( zztmp > zrho3 ) zrho10_3(ji,jj) = zzdep ! > 0.03 237 IF( zztmp > zdelr(ji,jj) ) zpycn (ji,jj) = zzdep ! > equi. delta T(10m) - 0.2 238 ! 239 END DO 240 END DO 241 END DO 242 243 CALL iom_put( "mld_dt02", zabs2 ) ! MLD abs(delta t) - 0.2 244 CALL iom_put( "topthdep", ztm2 ) ! T(10) - 0.2 245 CALL iom_put( "mldr10_3", zrho10_3 ) ! MLD delta rho(10m) = 0.03 246 CALL iom_put( "pycndep" , zpycn ) ! MLD delta rho equi. delta T(10m) = 0.2 247 CALL iom_put( "tinv" , ztinv ) ! max. temp. inv. (t10 ref) 248 CALL iom_put( "depti" , zdepinv ) ! depth of max. temp. inv. (t10 ref) 249 250 251 ! ----------------------------------- ! 252 ! search deepest level above 20C/28C ! 253 ! ----------------------------------- ! 254 ik20(:,:) = 1 255 ik28(:,:) = 1 256 DO jk = 1, jpkm1 ! beware temperature is not always decreasing with depth => loop from top to bottom 257 DO jj = 1, jpj 258 DO ji = 1, jpi 259 zztmp = tsn(ji,jj,jk,jp_tem) 260 IF( zztmp >= 20. ) ik20(ji,jj) = jk 261 IF( zztmp >= 28. ) ik28(ji,jj) = jk 262 END DO 263 END DO 264 END DO 265 266 ! --------------------------- ! 267 ! Depth of 20C/28C isotherm ! 268 ! --------------------------- ! 269 DO jj = 1, jpj 270 DO ji = 1, jpi 271 ! 272 zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1) ! depth of the oean bottom 273 ! 274 iid = ik20(ji,jj) 275 IF( iid /= 1 ) THEN 276 zztmp = gdept_n(ji,jj,iid ) & ! linear interpolation 277 & + ( gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid) ) & 278 & * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem) ) & 279 & / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 280 hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1) ! bound by the ocean depth 281 ELSE 282 hd20(ji,jj) = 0._wp 283 ENDIF 284 ! 285 iid = ik28(ji,jj) 286 IF( iid /= 1 ) THEN 287 zztmp = gdept_n(ji,jj,iid ) & ! linear interpolation 288 & + ( gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid) ) & 289 & * ( 28.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem) ) & 290 & / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 291 hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1) ! bound by the ocean depth 292 ELSE 293 hd28(ji,jj) = 0._wp 294 ENDIF 295 296 END DO 297 END DO 298 CALL iom_put( "20d", hd20 ) ! depth of the 20 isotherm 299 CALL iom_put( "28d", hd28 ) ! depth of the 28 isotherm 300 301 ! ----------------------------- ! 302 ! Heat content of first 300 m ! 303 ! ----------------------------- ! 304 305 ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...) 306 ilevel = 0 307 zthick_0 = 0._wp 308 DO jk = 1, jpkm1 309 zthick_0 = zthick_0 + e3t_1d(jk) 310 IF( zthick_0 < 300. ) ilevel = jk 311 END DO 312 ! surface boundary condition 313 IF( ln_linssh ) THEN ; zthick(:,:) = sshn(:,:) ; htc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 314 ELSE ; zthick(:,:) = 0._wp ; htc3(:,:) = 0._wp 315 ENDIF 316 ! integration down to ilevel 317 DO jk = 1, ilevel 318 zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk) 319 htc3 (:,:) = htc3 (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk) 320 END DO 321 ! deepest layer 322 zthick(:,:) = 300. - zthick(:,:) ! remaining thickness to reach 300m 323 DO jj = 1, jpj 324 DO ji = 1, jpi 325 htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) & 326 & * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1) 327 END DO 328 END DO 329 ! from temperature to heat contain 330 zcoef = rau0 * rcp 331 htc3(:,:) = zcoef * htc3(:,:) 332 CALL iom_put( "hc300", htc3 ) ! first 300m heat content 333 ! 334 IF( ln_timing ) CALL timing_stop('dia_hth') 335 ! 336 END SUBROUTINE dia_hth 337 338 #else 339 !!---------------------------------------------------------------------- 340 !! Default option : Empty module 341 !!---------------------------------------------------------------------- 342 LOGICAL , PUBLIC, PARAMETER :: lk_diahth = .FALSE. !: thermocline-20d depths flag 343 CONTAINS 344 SUBROUTINE dia_hth( kt ) ! Empty routine 345 IMPLICIT NONE 346 INTEGER, INTENT( in ) :: kt 347 WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt 348 END SUBROUTINE dia_hth 349 #endif 350 351 !!====================================================================== 334 END SUBROUTINE dia_hth_init 352 335 END MODULE diahth -
NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/step.F90
r11715 r12517 101 101 IF( kstp == nit000 ) THEN ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS) 102 102 CALL iom_init( cxios_context ) ! for model grid (including passible AGRIF zoom) 103 IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" ) ! for coarse grid 103 IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" ) ! for coarse grid 104 CALL dia_hth_init ! extra ML depth diagnostics, thermocline depths 104 105 ENDIF 105 106 IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init) … … 205 206 IF( ln_floats ) CALL flo_stp ( kstp ) ! drifting Floats 206 207 IF( ln_diacfl ) CALL dia_cfl ( kstp ) ! Courant number diagnostics 207 IF( l k_diahth ) CALL dia_hth ( kstp ) ! Thermocline depth (20 degres isotherm depth)208 IF( ll_diahth ) CALL dia_hth ( kstp ) ! Thermocline depth (20 degres isotherm depth) 208 209 IF( ln_diadct ) CALL dia_dct ( kstp ) ! Transports 209 210 CALL dia_ar5 ( kstp ) ! ar5 diag
Note: See TracChangeset
for help on using the changeset viewer.
