Changeset 15095
- Timestamp:
- 2021-07-07T13:22:23+02:00 (3 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_r4.2RC_GO8_package/src
- Files:
-
- 4 edited
-
OCE/LDF/ldftra.F90 (modified) (8 diffs)
-
OCE/SBC/sbcssm.F90 (modified) (2 diffs)
-
OCE/ZDF/zdfmxl.F90 (modified) (4 diffs)
-
SAS/sbcssm.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_r4.2RC_GO8_package/src/OCE/LDF/ldftra.F90
r14834 r15095 72 72 REAL(wp), PUBLIC :: rn_Ue !: lateral diffusive velocity [m/s] 73 73 REAL(wp), PUBLIC :: rn_Le !: lateral diffusive length [m] 74 INTEGER, PUBLIC :: nn_ldfeiv_shape !: shape of bounding coefficient (Treguier et al formulation only) 75 74 76 75 77 ! ! Flag to control the type of lateral diffusive operator … … 405 407 !!---------------------------------------------------------------------- 406 408 ! 407 IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21) THEN ! eddy induced velocity coefficients409 IF( ln_ldfeiv .AND. ( nn_aei_ijk_t == 21 .OR. nn_aei_ijk_t == 22 ) ) THEN ! eddy induced velocity coefficients 408 410 ! ! =F(growth rate of baroclinic instability) 409 411 ! ! max value aeiv_0 ; decreased to 0 within 20N-20S … … 496 498 !! 497 499 NAMELIST/namtra_eiv/ ln_ldfeiv , ln_ldfeiv_dia, & ! eddy induced velocity (eiv) 498 & nn_aei_ijk_t, rn_Ue, rn_Le ! eiv coefficient 500 & nn_aei_ijk_t, rn_Ue, rn_Le, & ! eiv coefficient 501 & nn_ldfeiv_shape 499 502 !!---------------------------------------------------------------------- 500 503 ! … … 517 520 WRITE(numout,*) ' eiv streamfunction & velocity diag. ln_ldfeiv_dia = ', ln_ldfeiv_dia 518 521 WRITE(numout,*) ' coefficients :' 519 WRITE(numout,*) ' type of time-space variation nn_aei_ijk_t = ', nn_a ht_ijk_t522 WRITE(numout,*) ' type of time-space variation nn_aei_ijk_t = ', nn_aei_ijk_t 520 523 WRITE(numout,*) ' lateral diffusive velocity (if cst) rn_Ue = ', rn_Ue, ' m/s' 521 524 WRITE(numout,*) ' lateral diffusive length (if cst) rn_Le = ', rn_Le, ' m' … … 589 592 IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )' 590 593 IF(lwp) WRITE(numout,*) ' maximum allowed value: aei0 = ', aei0, ' m2/s' 594 IF(lwp) WRITE(numout,*) ' shape of bounding coefficient : ',nn_ldfeiv_shape 591 595 ! 592 596 l_ldfeiv_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90 … … 637 641 ! 638 642 INTEGER :: ji, jj, jk ! dummy loop indices 639 REAL(wp) :: zfw, ze3w, zn2, z1_f20, zzaei ! local scalars 640 REAL(wp), DIMENSION(jpi,jpj) :: zn, zah, zhw, zRo, zaeiw ! 2D workspace 643 REAL(wp) :: zfw, ze3w, zn2, z1_f20, zzaei, z2_3 ! local scalars 644 REAL(wp), DIMENSION(jpi,jpj) :: zn, zah, zhw, zRo, zRo_lim, zTclinic_recip, zaeiw, zratio ! 2D workspace 645 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmodslp ! 3D workspace 641 646 !!---------------------------------------------------------------------- 642 647 ! 643 648 zn (:,:) = 0._wp ! Local initialization 649 zmodslp(:,:,:) = 0._wp 644 650 zhw(:,:) = 5._wp 645 651 zah(:,:) = 0._wp 646 652 zRo(:,:) = 0._wp 653 zRo_lim(:,:) = 0._wp 654 zTclinic_recip(:,:) = 0._wp 655 zratio(:,:) = 0._wp 656 zaeiw(:,:) = 0._wp 657 647 658 ! ! Compute lateral diffusive coefficient at T-point 648 659 IF( ln_traldf_triad ) THEN … … 671 682 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 672 683 ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 673 zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 674 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w 684 zmodslp(ji,jj,jk) = wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 685 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 686 zah(ji,jj) = zah(ji,jj) + zn2 * zmodslp(ji,jj,jk) * ze3w 675 687 zhw(ji,jj) = zhw(ji,jj) + ze3w 676 688 END_3D … … 680 692 zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 ) 681 693 ! Rossby radius at w-point taken betwenn 2 km and 40km 682 zRo(ji,jj) = MAX( 2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 ) ) 694 zRo(ji,jj) = .4 * zn(ji,jj) / zfw 695 zRo_lim(ji,jj) = MAX( 2.e3 , MIN( zRo(ji,jj), 40.e3 ) ) 683 696 ! Compute aeiw by multiplying Ro^2 and T^-1 684 zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1) 697 zTclinic_recip(ji,jj) = SQRT( MAX(zah(ji,jj),0._wp) / zhw(ji,jj) ) * tmask(ji,jj,1) 698 zaeiw(ji,jj) = zRo_lim(ji,jj) * zRo_lim(ji,jj) * zTclinic_recip(ji,jj) 685 699 END_2D 686 700 IF( iom_use('N_2d') ) CALL iom_put('N_2d',zn(:,:)/zhw(:,:)) 701 IF( iom_use('modslp') ) CALL iom_put('modslp',SQRT(zmodslp(:,:,:)) ) 702 CALL iom_put('RossRad',zRo) 703 CALL iom_put('RossRadlim',zRo_lim) 704 CALL iom_put('Tclinic_recip',zTclinic_recip) 687 705 ! !== Bound on eiv coeff. ==! 688 706 z1_f20 = 1._wp / ( 2._wp * omega * sin( rad * 20._wp ) ) 689 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 690 zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease 691 zaeiw(ji,jj) = MIN( zzaei , paei0 ) ! Max value = paei0 692 END_2D 707 z2_3 = 2._wp/3._wp 708 709 SELECT CASE(nn_ldfeiv_shape) 710 CASE(0) !! Standard shape applied - decrease in tropics and cap. 711 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 712 zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease 713 zaeiw(ji,jj) = MIN( zzaei, paei0 ) 714 END_2D 715 716 CASE(1) !! Abrupt cut-off on Rossby radius: 717 ! JD : modifications here to introduce scaling by local rossby radius of deformation vs local grid scale 718 ! arbitrary decision that GM is de-activated if local rossy radius larger than 2 times local grid scale 719 ! based on Hallberg (2013) 720 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 721 IF ( zRo(ji,jj) >= ( 2._wp * MIN( e1t(ji,jj), e2t(ji,jj) ) ) ) THEN 722 ! TODO : use a version of zRo that integrates over a few time steps ? 723 zaeiw(ji,jj) = 0._wp 724 ELSE 725 zaeiw(ji,jj) = MIN( zaeiw(ji,jj), paei0 ) 726 ENDIF 727 END_2D 728 CASE(2) !! Rossby radius ramp type 1: 729 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 730 zratio(ji,jj) = zRo(ji,jj)/MIN(e1t(ji,jj),e2t(ji,jj)) 731 zaeiw(ji,jj) = MIN( zaeiw(ji,jj), MAX( 0._wp, MIN( 1._wp, z2_3*(2._wp - zratio(ji,jj)) ) ) * paei0 ) 732 END_2D 733 CALL iom_put('RR_GS',zratio) 734 CASE(3) !! Rossby radius ramp type 2: 735 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 736 zratio(ji,jj) = MIN(e1t(ji,jj),e2t(ji,jj))/zRo(ji,jj) 737 zaeiw(ji,jj) = MIN( zaeiw(ji,jj), MAX( 0._wp, MIN( 1._wp, z2_3*( zratio(ji,jj) - 0.5_wp ) ) ) * paei0 ) 738 END_2D 739 CASE(4) !! bathymetry ramp: 740 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 741 zaeiw(ji,jj) = MIN( zaeiw(ji,jj), MAX( 0._wp, MIN( 1._wp, 0.001*(ht_0(ji,jj) - 2000._wp) ) ) * paei0 ) 742 END_2D 743 CASE(5) !! Rossby radius ramp type 1 applied to Treguier et al coefficient rather than cap: 744 !! Note the ramp is RR/GS=[2.0,1.0] (not [2.0,0.5] as for cases 2,3) and we ramp up 745 !! to 5% of the Treguier et al coefficient, aiming for peak values of around 100m2/s 746 !! at high latitudes rather than 2000m2/s which is what you get in eORCA025 with an 747 !! uncapped coefficient. 748 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 749 zratio(ji,jj) = zRo(ji,jj)/MIN(e1t(ji,jj),e2t(ji,jj)) 750 zaeiw(ji,jj) = MAX( 0._wp, MIN( 1._wp, 2._wp - zratio(ji,jj) ) ) * 0.05 * zaeiw(ji,jj) 751 zaeiw(ji,jj) = MIN( zaeiw(ji,jj), paei0 ) 752 END_2D 753 CALL iom_put('RR_GS',zratio) 754 CASE DEFAULT 755 CALL ctl_stop('ldf_eiv: Unrecognised option for nn_ldfeiv_shape.') 756 END SELECT 693 757 CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1.0_wp ) ! lateral boundary condition 694 758 ! -
NEMO/branches/UKMO/NEMO_r4.2RC_GO8_package/src/OCE/SBC/sbcssm.F90
r15062 r15095 59 59 REAL(wp) :: zcoef, zf_sbc ! local scalar 60 60 REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts 61 !!--------------------------------------------------------------------- 61 CHARACTER(len=4),SAVE :: stype 62 !!--------------------------------------------------------------------- 63 IF( kt == nit000 ) THEN 64 IF( ln_TEOS10 ) THEN 65 stype='abs' ! teos-10: using absolute salinity (sst is converted to potential temperature for the surface module) 66 ELSE IF( ln_EOS80 ) THEN 67 stype='pra' ! eos-80: using practical salinity 68 ELSE IF ( ln_SEOS) THEN 69 stype='seos' ! seos using Simplified Equation of state (sst is converted to potential temperature for the surface module) 70 ENDIF 71 ENDIF 62 72 ! 63 73 ! !* surface T-, U-, V- ocean level variables (T, S, depth, velocity) … … 170 180 CALL iom_put( 'ssu_m', ssu_m ) 171 181 CALL iom_put( 'ssv_m', ssv_m ) 172 CALL iom_put( 'sst_m ', sst_m )173 CALL iom_put( 'sss_m ', sss_m )182 CALL iom_put( 'sst_m_pot', sst_m ) 183 CALL iom_put( 'sss_m_'//stype, sss_m ) 174 184 CALL iom_put( 'ssh_m', ssh_m ) 175 185 CALL iom_put( 'e3t_m', e3t_m ) -
NEMO/branches/UKMO/NEMO_r4.2RC_GO8_package/src/OCE/ZDF/zdfmxl.F90
r14834 r15095 16 16 USE trc_oce , ONLY: l_offline ! ocean space and time domain variables 17 17 USE zdf_oce ! ocean vertical physics 18 USE eosbn2 ! for zdf_mxl_zint 18 19 ! 19 20 USE in_out_manager ! I/O manager … … 32 33 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlp !: mixed layer depth (rho=rho0+zdcrit) [m] (used by LDF) 33 34 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlpt !: depth of the last T-point inside the mixed layer [m] (used by LDF) 35 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hmld_zint !: vertically-interpolated mixed layer depth [m] 36 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: htc_mld ! Heat content of hmld_zint 37 LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ll_found ! Is T_b to be found by interpolation ? 38 LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ll_belowml ! Flag points below mixed layer when ll_found=F 34 39 35 40 REAL(wp), PUBLIC :: rho_c = 0.01_wp !: density criterion for mixed layer depth 36 41 REAL(wp), PUBLIC :: avt_c = 5.e-4_wp ! Kz criterion for the turbocline depth 42 43 TYPE, PUBLIC :: MXL_ZINT !: Structure for MLD defs 44 INTEGER :: mld_type ! mixed layer type 45 REAL(wp) :: zref ! depth of initial T_ref 46 REAL(wp) :: dT_crit ! Critical temp diff 47 REAL(wp) :: iso_frac ! Fraction of rn_dT_crit 48 END TYPE MXL_ZINT 37 49 38 50 !! * Substitutions … … 52 64 zdf_mxl_alloc = 0 ! set to zero if no array to be allocated 53 65 IF( .NOT. ALLOCATED( nmln ) ) THEN 54 ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc ) 66 ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj), & 67 & htc_mld(jpi,jpj), ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc ) 55 68 ! 56 69 CALL mpp_sum ( 'zdfmxl', zdf_mxl_alloc ) … … 117 130 ! 118 131 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ' ) 132 ! Vertically-interpolated mixed-layer depth diagnostic 133 CALL zdf_mxl_zint( kt , Kmm) 134 ! 119 135 ! 120 136 END SUBROUTINE zdf_mxl 137 138 SUBROUTINE zdf_mxl_zint_mld( sf , Kmm) 139 !!---------------------------------------------------------------------------------- 140 !! *** ROUTINE zdf_mxl_zint_mld *** 141 ! 142 ! Calculate vertically-interpolated mixed layer depth diagnostic. 143 ! 144 ! This routine can calculate the mixed layer depth diagnostic suggested by 145 ! Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate 146 ! vertically-interpolated mixed-layer depth diagnostics with other parameter 147 ! settings set in the namzdf_mldzint namelist. 148 ! 149 ! If mld_type=1 the mixed layer depth is calculated as the depth at which the 150 ! density has increased by an amount equivalent to a temperature difference of 151 ! 0.8C at the surface. 152 ! 153 ! For other values of mld_type the mixed layer is calculated as the depth at 154 ! which the temperature differs by 0.8C from the surface temperature. 155 ! 156 ! David Acreman, Daley Calvert 157 ! 158 !!----------------------------------------------------------------------------------- 159 160 TYPE(MXL_ZINT), INTENT(in) :: sf 161 INTEGER, INTENT(in) :: Kmm ! ocean time level index 162 163 ! Diagnostic criteria 164 INTEGER :: nn_mld_type ! mixed layer type 165 REAL(wp) :: rn_zref ! depth of initial T_ref 166 REAL(wp) :: rn_dT_crit ! Critical temp diff 167 REAL(wp) :: rn_iso_frac ! Fraction of rn_dT_crit used 168 169 ! Local variables 170 REAL(wp), PARAMETER :: zepsilon = 1.e-30 ! local small value 171 INTEGER, DIMENSION(jpi,jpj) :: ikmt ! number of active tracer levels 172 INTEGER, DIMENSION(jpi,jpj) :: ik_ref ! index of reference level 173 INTEGER, DIMENSION(jpi,jpj) :: ik_iso ! index of last uniform temp level 174 REAL, DIMENSION(jpi,jpj,jpk) :: zT ! Temperature or density 175 REAL, DIMENSION(jpi,jpj) :: ppzdep ! depth for use in calculating d(rho) 176 REAL, DIMENSION(jpi,jpj) :: zT_ref ! reference temperature 177 REAL :: zT_b ! base temperature 178 REAL, DIMENSION(jpi,jpj,jpk) :: zdTdz ! gradient of zT 179 REAL, DIMENSION(jpi,jpj,jpk) :: zmoddT ! Absolute temperature difference 180 REAL :: zdz ! depth difference 181 REAL :: zdT ! temperature difference 182 REAL, DIMENSION(jpi,jpj) :: zdelta_T ! difference critereon 183 REAL, DIMENSION(jpi,jpj) :: zRHO1, zRHO2 ! Densities 184 INTEGER :: ji, jj, jk ! loop counter 185 186 !!------------------------------------------------------------------------------------- 187 ! 188 ! Unpack structure 189 nn_mld_type = sf%mld_type 190 rn_zref = sf%zref 191 rn_dT_crit = sf%dT_crit 192 rn_iso_frac = sf%iso_frac 193 194 ! Set the mixed layer depth criterion at each grid point 195 IF( nn_mld_type == 0 ) THEN 196 zdelta_T(:,:) = rn_dT_crit 197 zT(:,:,:) = rhop(:,:,:) 198 ELSE IF( nn_mld_type == 1 ) THEN 199 ppzdep(:,:)=0.0 200 call eos ( ts(:,:,1,:,Kmm), ppzdep(:,:), zRHO1(:,:) ) 201 ! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST 202 ! [assumes number of tracers less than number of vertical levels] 203 zT(:,:,1:jpts)=ts(:,:,1,1:jpts,Kmm) 204 zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit 205 CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) ) 206 zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rho0 207 ! RHO from eos (2d version) doesn't calculate north or east halo: 208 CALL lbc_lnk( 'zdfmxl', zdelta_T, 'T', 1. ) 209 zT(:,:,:) = rhop(:,:,:) 210 ELSE 211 zdelta_T(:,:) = rn_dT_crit 212 zT(:,:,:) = ts(:,:,:,jp_tem,Kmm) 213 END IF 214 215 ! Calculate the gradient of zT and absolute difference for use later 216 DO jk = 1 ,jpk-2 217 zdTdz(:,:,jk) = ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w(:,:,jk+1,Kmm) 218 zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 219 END DO 220 221 ! Find density/temperature at the reference level (Kara et al use 10m). 222 ! ik_ref is the index of the box centre immediately above or at the reference level 223 ! Find rn_zref in the array of model level depths and find the ref 224 ! density/temperature by linear interpolation. 225 DO jk = jpkm1, 2, -1 226 WHERE ( gdept(:,:,jk,Kmm) > rn_zref ) 227 ik_ref(:,:) = jk - 1 228 zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - gdept(:,:,jk-1,Kmm) ) 229 END WHERE 230 END DO 231 232 ! If the first grid box centre is below the reference level then use the 233 ! top model level to get zT_ref 234 WHERE ( gdept(:,:,1,Kmm) > rn_zref ) 235 zT_ref = zT(:,:,1) 236 ik_ref = 1 237 END WHERE 238 239 ! The number of active tracer levels is 1 less than the number of active w levels 240 ikmt(:,:) = mbkt(:,:) - 1 241 242 ! Initialize / reset 243 ll_found(:,:) = .false. 244 245 IF ( rn_iso_frac - zepsilon > 0. ) THEN 246 ! Search for a uniform density/temperature region where adjacent levels 247 ! differ by less than rn_iso_frac * deltaT. 248 ! ik_iso is the index of the last level in the uniform layer 249 ! ll_found indicates whether the mixed layer depth can be found by interpolation 250 ik_iso(:,:) = ik_ref(:,:) 251 DO jj = 1, jpj ! Changed from nlcj 252 DO ji = 1, jpi ! Changed from nlci 253 !CDIR NOVECTOR 254 DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1 255 IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN 256 ik_iso(ji,jj) = jk 257 ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 258 EXIT 259 END IF 260 END DO 261 END DO 262 END DO 263 264 ! Use linear interpolation to find depth of mixed layer base where possible 265 hmld_zint(:,:) = rn_zref 266 DO jj = 1, jpj 267 DO ji = 1, jpi 268 IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN 269 zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 270 hmld_zint(ji,jj) = gdept(ji,jj,ik_iso(ji,jj),Kmm) + zdz 271 END IF 272 END DO 273 END DO 274 END IF 275 276 ! If ll_found = .false. then calculate MLD using difference of zdelta_T 277 ! from the reference density/temperature 278 279 ! Prevent this section from working on land points 280 WHERE ( tmask(:,:,1) /= 1.0 ) 281 ll_found = .true. 282 END WHERE 283 284 DO jk=1, jpk 285 ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:) 286 END DO 287 288 ! Set default value where interpolation cannot be used (ll_found=false) 289 DO jj = 1, jpj 290 DO ji = 1, jpi 291 IF ( .not. ll_found(ji,jj) ) hmld_zint(ji,jj) = gdept(ji,jj,ikmt(ji,jj),Kmm) 292 END DO 293 END DO 294 295 DO jj = 1, jpj 296 DO ji = 1, jpi 297 !CDIR NOVECTOR 298 DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj) 299 IF ( ll_found(ji,jj) ) EXIT 300 IF ( ll_belowml(ji,jj,jk) ) THEN 301 zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) ) 302 zdT = zT_b - zT(ji,jj,jk-1) 303 zdz = zdT / zdTdz(ji,jj,jk-1) 304 hmld_zint(ji,jj) = gdept(ji,jj,jk-1,Kmm) + zdz 305 EXIT 306 END IF 307 END DO 308 END DO 309 END DO 310 311 hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1) 312 ! 313 END SUBROUTINE zdf_mxl_zint_mld 314 315 SUBROUTINE zdf_mxl_zint_htc( kt , Kmm) 316 !!---------------------------------------------------------------------- 317 !! *** ROUTINE zdf_mxl_zint_htc *** 318 !! 319 !! ** Purpose : 320 !! 321 !! ** Method : 322 !!---------------------------------------------------------------------- 323 324 INTEGER, INTENT(in) :: kt ! ocean time-step index 325 INTEGER, INTENT(in) :: Kmm ! ocean time level index 326 327 INTEGER :: ji, jj, jk 328 INTEGER :: ikmax 329 REAL(wp) :: zc, zcoef 330 ! 331 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilevel 332 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zthick_0, zthick 333 334 !!---------------------------------------------------------------------- 335 336 IF( .NOT. ALLOCATED(ilevel) ) THEN 337 ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), & 338 & zthick(jpi,jpj), STAT=ji ) 339 IF( lk_mpp ) CALL mpp_sum( 'zdfmxl', ji ) 340 IF( ji /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' ) 341 ENDIF 342 343 ! Find last whole model T level above the MLD 344 ilevel(:,:) = 0 345 zthick_0(:,:) = 0._wp 346 347 DO jk = 1, jpkm1 348 DO jj = 1, jpj 349 DO ji = 1, jpi 350 zthick_0(ji,jj) = zthick_0(ji,jj) + e3t(ji,jj,jk,Kmm) 351 IF( zthick_0(ji,jj) < hmld_zint(ji,jj) ) ilevel(ji,jj) = jk 352 END DO 353 END DO 354 WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2) 355 WRITE(numout,*) 'gdepw(jk+1 =',jk+1,') =',gdepw(2,2,jk+1,Kmm) 356 END DO 357 358 ! Surface boundary condition 359 IF( ln_linssh ) THEN ; zthick(:,:) = ssh(:,:,Kmm) ; htc_mld(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) * tmask(:,:,1) 360 ELSE ; zthick(:,:) = 0._wp ; htc_mld(:,:) = 0._wp 361 ENDIF 362 363 ! Deepest whole T level above the MLD 364 ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 ) 365 366 ! Integration down to last whole model T level 367 DO jk = 1, ikmax 368 DO jj = 1, jpj 369 DO ji = 1, jpi 370 zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1 ) ) ! 0 below ilevel 371 zthick(ji,jj) = zthick(ji,jj) + zc 372 htc_mld(ji,jj) = htc_mld(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk) 373 END DO 374 END DO 375 END DO 376 377 ! Subsequent partial T level 378 zthick(:,:) = hmld_zint(:,:) - zthick(:,:) ! remaining thickness to reach MLD 379 380 DO jj = 1, jpj 381 DO ji = 1, jpi 382 htc_mld(ji,jj) = htc_mld(ji,jj) + ts(ji,jj,ilevel(ji,jj)+1,jp_tem,Kmm) & 383 & * MIN( e3t(ji,jj,ilevel(ji,jj)+1,Kmm), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1) 384 END DO 385 END DO 386 387 WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2) 388 389 ! Convert to heat content 390 zcoef = rho0 * rcp 391 htc_mld(:,:) = zcoef * htc_mld(:,:) 392 393 END SUBROUTINE zdf_mxl_zint_htc 394 395 SUBROUTINE zdf_mxl_zint( kt , Kmm) 396 !!---------------------------------------------------------------------- 397 !! *** ROUTINE zdf_mxl_zint *** 398 !! 399 !! ** Purpose : 400 !! 401 !! ** Method : 402 !!---------------------------------------------------------------------- 403 404 INTEGER, INTENT(in) :: kt ! ocean time-step index 405 INTEGER, INTENT(in) :: Kmm ! ocean time level index 406 407 INTEGER :: ios 408 INTEGER :: jn 409 410 INTEGER :: nn_mld_diag = 0 ! number of diagnostics 411 412 CHARACTER(len=1) :: cmld 413 414 TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 415 TYPE(MXL_ZINT), SAVE, DIMENSION(5) :: mld_diags 416 417 NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 418 419 !!---------------------------------------------------------------------- 420 421 IF( kt == nit000 ) THEN 422 READ ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901) 423 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist' ) 424 425 READ ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 ) 426 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist' ) 427 IF(lwm) WRITE ( numond, namzdf_mldzint ) 428 429 IF( nn_mld_diag > 5 ) CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' ) 430 431 mld_diags(1) = sn_mld1 432 mld_diags(2) = sn_mld2 433 mld_diags(3) = sn_mld3 434 mld_diags(4) = sn_mld4 435 mld_diags(5) = sn_mld5 436 437 IF( lwp .AND. (nn_mld_diag > 0) ) THEN 438 WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================' 439 WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)' 440 DO jn = 1, nn_mld_diag 441 WRITE(numout,*) 'MLD criterion',jn,':' 442 WRITE(numout,*) ' nn_mld_type =', mld_diags(jn)%mld_type 443 WRITE(numout,*) ' rn_zref =' , mld_diags(jn)%zref 444 WRITE(numout,*) ' rn_dT_crit =' , mld_diags(jn)%dT_crit 445 WRITE(numout,*) ' rn_iso_frac =', mld_diags(jn)%iso_frac 446 END DO 447 WRITE(numout,*) '====================================================================' 448 ENDIF 449 ENDIF 450 451 IF( nn_mld_diag > 0 ) THEN 452 DO jn = 1, nn_mld_diag 453 WRITE(cmld,'(I1)') jn 454 IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN 455 CALL zdf_mxl_zint_mld( mld_diags(jn) , Kmm) 456 457 IF( iom_use( "mldzint_"//cmld ) ) THEN 458 CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) ) 459 ENDIF 460 461 IF( iom_use( "mldhtc_"//cmld ) ) THEN 462 CALL zdf_mxl_zint_htc( kt , Kmm ) 463 CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:) ) 464 ENDIF 465 ENDIF 466 END DO 467 ENDIF 468 469 END SUBROUTINE zdf_mxl_zint 121 470 122 471 -
NEMO/branches/UKMO/NEMO_r4.2RC_GO8_package/src/SAS/sbcssm.F90
r15023 r15095 80 80 REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation 81 81 REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation 82 CHARACTER(len=4),SAVE :: stype 82 83 !!---------------------------------------------------------------------- 83 84 ! 84 85 IF( ln_timing ) CALL timing_start( 'sbc_ssm') 86 IF( kt == nit000 ) THEN 87 IF( ln_TEOS10 ) THEN 88 stype='abs' ! teos-10: using absolute salinity (sst is converted to potential temperature for the surface module) 89 ELSE IF( ln_EOS80 ) THEN 90 stype='pra' ! eos-80: using practical salinity 91 ELSE IF ( ln_SEOS) THEN 92 stype='seos' ! seos using Simplified Equation of state (sst is converted to potential temperature for the surface module) 93 ENDIF 94 ENDIF 85 95 86 96 IF ( l_sasread ) THEN … … 157 167 CALL iom_put( 'ssu_m', ssu_m ) 158 168 CALL iom_put( 'ssv_m', ssv_m ) 159 CALL iom_put( 'sst_m ', sst_m )160 CALL iom_put( 'sss_m ', sss_m )169 CALL iom_put( 'sst_m_pot', sst_m ) 170 CALL iom_put( 'sss_m_'//stype, sss_m ) 161 171 CALL iom_put( 'ssh_m', ssh_m ) 162 172 IF( .NOT.ln_linssh ) CALL iom_put( 'e3t_m', e3t_m )
Note: See TracChangeset
for help on using the changeset viewer.
