Changeset 15603 for branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
- Timestamp:
- 2021-12-16T10:39:55+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r12555 r15603 51 51 USE par_ice_2 52 52 #endif 53 USE stopack 53 54 54 55 IMPLICIT NONE … … 89 90 REAL(wp) :: rn_efac ! multiplication factor for evaporation (clem) 90 91 REAL(wp) :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem) 92 REAL(wp), ALLOCATABLE, SAVE :: rn_vfac0(:,:) ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem) 91 93 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 92 94 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 95 REAL(wp), PUBLIC :: rn_sfac ! multiplication factor for snow precipitation over sea-ice 93 96 94 97 !! * Substitutions … … 151 154 & sn_wndi, sn_wndj, sn_humi , sn_qsr , & 152 155 & sn_qlw , sn_tair, sn_prec , sn_snow, & 153 & sn_tdif, rn_zqt, rn_zu 156 & sn_tdif, rn_zqt, rn_zu, rn_sfac 154 157 !!--------------------------------------------------------------------- 155 158 ! … … 158 161 ! ! ====================== ! 159 162 ! 163 rn_sfac = 1._wp ! Default to one if missing from namelist 160 164 REWIND( numnam_ref ) ! Namelist namsbc_core in reference namelist : CORE bulk parameters 161 165 READ ( numnam_ref, namsbc_core, IOSTAT = ios, ERR = 901) … … 196 200 sfx(:,:) = 0._wp ! salt flux; zero unless ice is present (computed in limsbc(_2).F90) 197 201 ! 198 ENDIF 202 ALLOCATE ( rn_vfac0(jpi,jpj) ) 203 rn_vfac0(:,:) = rn_vfac 204 ! 205 ENDIF 206 207 #if defined key_traldf_c2d || key_traldf_c3d 208 IF( ln_stopack .AND. nn_spp_relw > 0 ) THEN 209 rn_vfac0(:,:) = rn_vfac 210 CALL spp_gen(kt, rn_vfac0, nn_spp_relw, rn_relw_sd, jk_spp_relw ) 211 ENDIF 212 #else 213 IF ( ln_stopack .AND. nn_spp_relw > 0 ) & 214 & CALL ctl_stop( 'sbc_blk_core: parameter perturbation will only work with '// & 215 'key_traldf_c2d or key_traldf_c3d') 216 #endif 199 217 200 218 CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step … … 205 223 #if defined key_cice 206 224 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 207 qlw_ice(:,:,1) = sf(jp_qlw)%fnow(:,:,1) 225 qlw_ice(:,:,1) = sf(jp_qlw)%fnow(:,:,1) 208 226 IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 209 227 ELSE ; qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 210 228 ENDIF 211 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 229 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 212 230 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 213 231 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac … … 219 237 ! 220 238 END SUBROUTINE sbc_blk_core 221 222 239 240 223 241 SUBROUTINE blk_oce_core( kt, sf, pst, pu, pv ) 224 242 !!--------------------------------------------------------------------- … … 230 248 !! ** Method : CORE bulk formulea for the ocean using atmospheric 231 249 !! fields read in sbc_read 232 !! 250 !! 233 251 !! ** Outputs : - utau : i-component of the stress at U-point (N/m2) 234 252 !! - vtau : j-component of the stress at V-point (N/m2) … … 268 286 ! local scalars ( place there for vector optimisation purposes) 269 287 zcoef_qsatw = 0.98 * 640380. / rhoa 270 288 271 289 zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 272 290 … … 276 294 277 295 ! ... components ( U10m - U_oce ) at T-point (unmasked) 278 zwnd_i(:,:) = 0.e0 296 zwnd_i(:,:) = 0.e0 279 297 zwnd_j(:,:) = 0.e0 280 298 #if defined key_cyclone … … 289 307 DO jj = 2, jpjm1 290 308 DO ji = fs_2, fs_jpim1 ! vect. opt. 291 zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) )292 zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) )309 zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) ) 310 zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) ) 293 311 END DO 294 312 END DO … … 320 338 CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm, & 321 339 & Cd, Ch, Ce, zt_zu, zq_zu ) 322 340 323 341 ! ... tau module, i and j component 324 342 DO jj = 1, jpj … … 351 369 CALL lbc_lnk( vtau(:,:), 'V', -1. ) 352 370 353 371 354 372 ! Turbulent fluxes over ocean 355 373 ! ----------------------------- … … 376 394 CALL prt_ctl( tab2d_1=zst , clinfo1=' blk_oce_core: zst : ') 377 395 ENDIF 378 396 379 397 ! ----------------------------------------------------------------------------- ! 380 398 ! III Total FLUXES ! … … 384 402 & - sf(jp_prec)%fnow(:,:,1) * rn_pfac ) * tmask(:,:,1) 385 403 ! 386 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar 404 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar 387 405 & - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus & ! remove latent melting heat for solid precip 388 406 & - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST … … 425 443 ! 426 444 END SUBROUTINE blk_oce_core 427 428 445 446 429 447 #if defined key_lim2 || defined key_lim3 430 448 SUBROUTINE blk_ice_core_tau … … 465 483 DO ji = 2, jpim1 ! B grid : NO vector opt 466 484 ! ... scalar wind at I-point (fld being at T-point) 467 zwndi_f = 0.25 * ( sf(jp_wndi)%fnow(ji-1,jj ,1) + sf(jp_wndi)%fnow(ji ,jj ,1) & 468 & + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji ,jj-1,1) ) - rn_vfac * u_ice(ji,jj) 469 zwndj_f = 0.25 * ( sf(jp_wndj)%fnow(ji-1,jj ,1) + sf(jp_wndj)%fnow(ji ,jj ,1) & 470 & + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji ,jj-1,1) ) - rn_vfac * v_ice(ji,jj) 485 zwndi_f = 0.25 * ( sf(jp_wndi)%fnow(ji-1,jj ,1) + sf(jp_wndi)%fnow(ji ,jj ,1) & 486 & + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji ,jj-1,1) ) & 487 & - rn_vfac0(ji,jj) * u_ice(ji,jj) 488 zwndj_f = 0.25 * ( sf(jp_wndj)%fnow(ji-1,jj ,1) + sf(jp_wndj)%fnow(ji ,jj ,1) & 489 & + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji ,jj-1,1) ) & 490 & - rn_vfac0(ji,jj) * v_ice(ji,jj) 471 491 zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 472 492 ! ... ice stress at I-point … … 474 494 vtau_ice(ji,jj) = zwnorm_f * zwndj_f 475 495 ! ... scalar wind at T-point (fld being at T-point) 476 zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( u_ice(ji,jj+1) + u_ice(ji+1,jj+1) & 477 & + u_ice(ji,jj ) + u_ice(ji+1,jj ) ) 478 zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( v_ice(ji,jj+1) + v_ice(ji+1,jj+1) & 479 & + v_ice(ji,jj ) + v_ice(ji+1,jj ) ) 496 zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) & 497 & - rn_vfac0(ji,jj) * 0.25 * ( u_ice(ji,jj+1) + u_ice(ji+1,jj+1) & 498 & + u_ice(ji,jj ) + u_ice(ji+1,jj ) ) 499 zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) & 500 & - rn_vfac0(ji,jj) * 0.25 * ( v_ice(ji,jj+1) + v_ice(ji+1,jj+1) & 501 & + v_ice(ji,jj ) + v_ice(ji+1,jj ) ) 480 502 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 481 503 END DO … … 488 510 DO jj = 2, jpj 489 511 DO ji = fs_2, jpi ! vect. opt. 490 zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) )491 zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) )512 zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) ) 513 zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) ) 492 514 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 493 515 END DO … … 495 517 DO jj = 2, jpjm1 496 518 DO ji = fs_2, fs_jpim1 ! vect. opt. 497 utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) & 498 & * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 499 vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) & 500 & * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 519 utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) & 520 & * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) & 521 & - rn_vfac0(ji,jj) * u_ice(ji,jj) ) 522 vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) & 523 & * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) & 524 & - rn_vfac0(ji,jj) * v_ice(ji,jj) ) 501 525 END DO 502 526 END DO … … 513 537 514 538 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_core_tau') 515 539 516 540 END SUBROUTINE blk_ice_core_tau 517 541 … … 526 550 !! between atmosphere and sea-ice using CORE bulk 527 551 !! formulea, ice variables and read atmmospheric fields. 528 !! 552 !! 529 553 !! caution : the net upward water flux has with mm/day unit 530 554 !!--------------------------------------------------------------------- … … 546 570 IF( nn_timing == 1 ) CALL timing_start('blk_ice_core_flx') 547 571 ! 548 CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 572 CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 549 573 550 574 ! local scalars ( place there for vector optimisation purposes) … … 569 593 z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 570 594 ! lw sensitivity 571 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 595 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 572 596 573 597 ! ----------------------------! … … 579 603 z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 580 604 ! Latent Heat 581 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls * Cice * wndm_ice(ji,jj) & 605 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls * Cice * wndm_ice(ji,jj) & 582 606 & * ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1) ) ) 583 607 ! Latent heat sensitivity for ice (Dqla/Dt) … … 610 634 611 635 #if defined key_lim3 612 CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 636 CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 613 637 614 638 ! --- evaporation --- ! … … 620 644 ! --- evaporation minus precipitation --- ! 621 645 zsnw(:,:) = 0._wp 622 CALL lim_thd_snwblow( pfrld, zsnw ) ! snow distribution over ice after wind blowing 646 CALL lim_thd_snwblow( pfrld, zsnw ) ! snow distribution over ice after wind blowing 623 647 emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 624 648 emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw … … 643 667 DO jl = 1, jpl 644 668 qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) ) 645 ! But we do not have Tice => consider it at 0 °C => evap=0669 ! But we do not have Tice => consider it at 0 degC => evap=0 646 670 END DO 647 671 648 CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 672 CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 649 673 #endif 650 674 … … 670 694 ! 671 695 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_core_flx') 672 696 673 697 END SUBROUTINE blk_ice_core_flx 674 698 #endif … … 683 707 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 684 708 !! 685 !! ** Method : Monin Obukhov Similarity Theory 709 !! ** Method : Monin Obukhov Similarity Theory 686 710 !! + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10) 687 711 !! … … 693 717 !! - better first guess of stability by checking air-sea difference of virtual temperature 694 718 !! rather than temperature difference only... 695 !! - added function "cd_neutral_10m" that uses the improved parametrization of 719 !! - added function "cd_neutral_10m" that uses the improved parametrization of 696 720 !! Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions! 697 721 !! - using code-wide physical constants defined into "phycst.mod" rather than redifining them … … 728 752 729 753 IF( nn_timing == 1 ) CALL timing_start('turb_core_2z') 730 754 731 755 CALL wrk_alloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd ) 732 756 CALL wrk_alloc( jpi,jpj, zeta_u, stab ) … … 740 764 U_zu = MAX( 0.5 , dU ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 741 765 742 !! First guess of stability: 766 !! First guess of stability: 743 767 ztmp0 = T_zt*(1. + 0.608*q_zt) - sst*(1. + 0.608*q_sat) ! air-sea difference of virtual pot. temp. at zt 744 768 stab = 0.5 + sign(0.5,ztmp0) ! stab = 1 if dTv > 0 => STABLE, 0 if unstable … … 754 778 Ce_n10 = 1.e-3*( 34.6 * sqrt_Cd_n10 ) 755 779 Ch_n10 = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) 756 780 757 781 !! Initializing transf. coeff. with their first guess neutral equivalents : 758 782 Cd = ztmp0 ; Ce = Ce_n10 ; Ch = Ch_n10 ; sqrt_Cd = sqrt_Cd_n10 … … 765 789 ! 766 790 ztmp1 = T_zu - sst ! Updating air/sea differences 767 ztmp2 = q_zu - q_sat 791 ztmp2 = q_zu - q_sat 768 792 769 793 ! Updating turbulent scales : (L&Y 2004 eq. (7)) 770 794 ztmp1 = Ch/sqrt_Cd*ztmp1 ! theta* 771 795 ztmp2 = Ce/sqrt_Cd*ztmp2 ! q* 772 796 773 797 ztmp0 = T_zu*(1. + 0.608*q_zu) ! virtual potential temperature at zu 774 798 775 799 ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 776 ztmp0 = (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu) 800 ztmp0 = (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu) 777 801 ! ( Cd*U_zu*U_zu is U*^2 at zu) 778 802 … … 781 805 zpsi_h_u = psi_h( zeta_u ) 782 806 zpsi_m_u = psi_m( zeta_u ) 783 807 784 808 !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c)) 785 809 IF ( .NOT. l_zt_equal_zu ) THEN … … 790 814 q_zu = max(0., q_zu) 791 815 END IF 792 816 793 817 IF( ln_cdgw ) THEN ! surface wave case 794 sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 818 sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 795 819 Cd = sqrt_Cd * sqrt_Cd 796 820 ELSE … … 802 826 ztmp0 = cd_neutral_10m(ztmp0) ! Cd_n10 803 827 sqrt_Cd_n10 = sqrt(ztmp0) 804 828 805 829 Ce_n10 = 1.e-3 * (34.6 * sqrt_Cd_n10) ! L&Y 2004 eq. (6b) 806 830 stab = 0.5 + sign(0.5,zeta_u) ! update stability … … 809 833 !! Update of transfer coefficients: 810 834 ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - zpsi_m_u) ! L&Y 2004 eq. (10a) 811 Cd = ztmp0 / ( ztmp1*ztmp1 ) 835 Cd = ztmp0 / ( ztmp1*ztmp1 ) 812 836 sqrt_Cd = SQRT( Cd ) 813 837 ENDIF … … 815 839 ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 816 840 ztmp2 = sqrt_Cd / sqrt_Cd_n10 817 ztmp1 = 1. + Ch_n10*ztmp0 841 ztmp1 = 1. + Ch_n10*ztmp0 818 842 Ch = Ch_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 819 843 ! 820 ztmp1 = 1. + Ce_n10*ztmp0 844 ztmp1 = 1. + Ce_n10*ztmp0 821 845 Ce = Ce_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 822 846 ! … … 836 860 FUNCTION cd_neutral_10m( zw10 ) 837 861 !!---------------------------------------------------------------------- 838 !! Estimate of the neutral drag coefficient at 10m as a function 862 !! Estimate of the neutral drag coefficient at 10m as a function 839 863 !! of neutral wind speed at 10m 840 864 !! … … 842 866 !! 843 867 !! Author: L. Brodeau, june 2014 844 !!---------------------------------------------------------------------- 868 !!---------------------------------------------------------------------- 845 869 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zw10 ! scalar wind speed at 10m (m/s) 846 870 REAL(wp), DIMENSION(jpi,jpj) :: cd_neutral_10m 847 871 ! 848 872 REAL(wp), DIMENSION(:,:), POINTER :: rgt33 849 !!---------------------------------------------------------------------- 873 !!---------------------------------------------------------------------- 850 874 ! 851 875 CALL wrk_alloc( jpi,jpj, rgt33 ) 852 876 ! 853 877 !! When wind speed > 33 m/s => Cyclone conditions => special treatment 854 rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) ) ! If zw10 < 33. => 0, else => 1 878 rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) ) ! If zw10 < 33. => 0, else => 1 855 879 cd_neutral_10m = 1.e-3 * ( & 856 880 & (1._wp - rgt33)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33.
Note: See TracChangeset
for help on using the changeset viewer.
