Changeset 606
- Timestamp:
- 2007-02-20T17:48:56+01:00 (17 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/SBC/flx_core.h90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/SBC/flx_core.h90
r472 r606 22 22 CHARACTER (LEN=32), DIMENSION (jpfile) :: clvarname 23 23 CHARACTER (LEN=50), DIMENSION (jpfile) :: clname 24 CHARACTER (LEN=32) :: clvarkatax , clvarkatay ! variable name for katabatic masks 25 CHARACTER (LEN=256) :: clnamekata ! filename for katabatic masks 26 24 27 INTEGER :: isnap 25 28 INTEGER, DIMENSION(jpfile) :: & … … 29 32 freqh ! Frequency for each forcing file (hours) 30 33 ! ! a negative value of -12 corresponds to monthly 34 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! used for modif wind stress in the first sea points 35 & rmskkatax, rmskkatay !: mask ocean to increase wind stress in first sea points 31 36 REAL(wp), DIMENSION(jpi,jpj,jpfile,2) :: & 32 37 flxdta !: set of NCAR 6hourly/daily/monthly fluxes 38 LOGICAL :: & 39 & ln_2m = .FALSE., & !: logical flag for height of air temp. and hum. 40 & ln_kata = .FALSE. !: logical flag for Katabatic winds enhancement. 33 41 !!---------------------------------------------------------------------- 34 42 !! History : … … 40 48 !! ! - Implement reading of 6-hourly fields 41 49 !! ! 06-02 (S. Masson, G. Madec) IOM read + NEMO "style" 50 !! ! 07-06 (P. Mathiot, J-M Molines) Katabatic winds enhancement 51 !! ! 12-06 (L. Brodeau) handle both 2m and 10m input fields 42 52 !!---------------------------------------------------------------------- 43 53 !! OPA 9.0 , LOCEAN-IPSL (2006) … … 111 121 INTEGER :: ihour ! current hour of the day at which we read the forcings 112 122 INTEGER :: & 113 imois, imois2, itime,& ! temporary integers114 i15 , iman , 123 imois, imois2, & ! temporary integers 124 i15 , iman , inum , & ! " " 115 125 ifpr , jfpr , & ! frequency of index for print in prihre 116 jj , ji! Loop indices126 jj , ji ! Loop indices 117 127 REAL(wp) :: & 118 128 zadatrjb, & ! date (noninteger day) at which we read the forcings 119 zsecond, & ! temporary scalars120 129 zxy , & ! scalar for temporal interpolation 121 130 zcof ! scalar 122 REAL(wp) :: em , aa123 REAL(wp) :: zind1, zind2, zind3124 131 REAL(wp), DIMENSION(jpi,jpj, jpfile) :: & 125 132 flxnow ! input flux values at current time step 126 133 REAL(wp), DIMENSION(jpi,jpj) :: & 127 sstk , t04128 REAL(wp), DIMENSION(jpi,jpj) :: &129 qtune , & ! artifical field for tuning q130 134 dqlw_ice , & ! long wave heat flx sensitivity over ice 131 135 dqsb_ice , & ! sensible heat flx sensitivity over ice 132 alb_ice , & ! albedo of ice133 136 alb_oce_os, & ! albedo of the ocean under overcast sky 134 137 alb_ice_os, & ! albedo of the ice under overcast sky 135 138 alb_ice_cs, & ! albedo of ice under clear sky 136 139 alb_oce_cs, & ! albedo of the ocean under clear sky 137 zsstlev, & ! SST from Levitus in K138 140 zsst, & ! nfbulk : mean SST 139 141 zsss, & ! nfbulk : mean tn_ice(:,:,1) … … 146 148 Ch, & ! coefficient for sensible heat flux 147 149 Ce, & ! coefficient for latent heat flux 148 Cd ! coefficient for wind stress 150 Cd, & ! coefficient for wind stress 151 zt_zu, & !: air temperature at wind speed height 152 zq_zu !: air spec. hum. at wind speed height 149 153 REAL(wp), PARAMETER :: & 150 154 rhoa = 1.22, & ! air density … … 157 161 catm1 ! 158 162 159 NAMELIST /namcore/ clname, clvarname, freqh163 NAMELIST /namcore/ clname, clvarname, freqh, ln_2m, ln_kata 160 164 !!--------------------------------------------------------------------- 161 165 … … 200 204 WRITE(numout,*)' flxmod/flx_core : global CORE fields in NetCDF format' 201 205 WRITE(numout,*)' ~~~~~~~~~~~~~~~ ' 206 WRITE(numout,*) ' ln_2m = ', ln_2m 207 WRITE(numout,*) ' ln_kata = ', ln_kata 202 208 WRITE(numout,*) ' list of files and frequency (hour), or monthly (-12) ' 203 209 DO ji = 1, jpfile … … 450 456 qsat(:,:) = 11637800*exp(-5897.8/zsss(:,:))/rhoa 451 457 ! 452 !--- NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point (NCAR Bulk formulae) 453 CALL TURB_CORE( 10., zsst(:,:), flxnow(:,:,7), qsatw(:,:), flxnow(:,:,4), & 454 & dUnormt(:,:), Cd(:,:), Ch(:,:), Ce(:,:) ) 455 458 459 ! CORE iterartive algo for computation of Cd, Ch, Ce at T-point : 460 ! =============================================================== 461 IF( ln_2m ) THEN 462 IF( kt == nit000 .AND. lwp ) THEN 463 WRITE(numout,*) 464 WRITE(numout,*)' flx_core : global CORE fields in NetCDF format' 465 WRITE(numout,*)' ~~~~~~~~~ ' 466 WRITE(numout,*) ' Calling TURB_CORE_2Z for bulk transfert coefficients' 467 WRITE(numout,*) 468 ENDIF 469 !! If air temp. and spec. hum. are given at different height (2m) than wind (10m) : 470 CALL TURB_CORE_2Z(2.,10.,zsst(:,:), flxnow(:,:,7), qsatw(:,:), flxnow(:,:,4), & 471 & dUnormt(:,:), Cd(:,:), Ch(:,:), Ce(:,:), zt_zu(:,:), zq_zu(:,:)) 472 ELSE 473 IF( kt == nit000 .AND. lwp ) THEN 474 WRITE(numout,*) 475 WRITE(numout,*)' flx_core : global CORE fields in NetCDF format' 476 WRITE(numout,*)' ~~~~~~~~~~ ' 477 WRITE(numout,*) ' Calling TURB_CORE_1Z for bulk transfert coefficients' 478 WRITE(numout,*) 479 ENDIF 480 !! If air temp. and spec. hum. are given at same height than wind (10m) : 481 CALL TURB_CORE_1Z(10., zsst(:,:), flxnow(:,:,7), qsatw(:,:), flxnow(:,:,4), & 482 & dUnormt(:,:), Cd(:,:), Ch(:,:), Ce(:,:) ) 483 ENDIF 484 456 485 ! II.1) Momentum over ocean and ice 457 486 ! --------------------------------- … … 466 495 CALL lbc_lnk( tauyt(:,:), 'T', -1. ) 467 496 ! 497 498 ! Katabatic winds parametrization 499 ! ------------------------------- 500 IF( ln_kata ) THEN 501 ! 502 IF( kt == nit000 ) THEN 503 IF (lwp ) WRITE(numout,*) ' Katabatic winds parametrization is active' 504 clnamekata = 'katamask.nc' 505 clvarkatax = 'katamaskx' 506 clvarkatay = 'katamasky' 507 508 #if defined key_agrif 509 IF ( .NOT. Agrif_Root() ) THEN 510 clnamekata = TRIM(Agrif_CFixed())//'_'//TRIM(clnamekata) 511 ENDIF 512 #endif 513 CALL iom_open( TRIM(clnamekata), inum ) 514 CALL iom_get ( inum, jpdom_data, TRIM(clvarkatax), rmskkatax ) 515 CALL iom_get ( inum, jpdom_data, TRIM(clvarkatay), rmskkatay ) 516 CALL iom_close ( inum ) 517 518 WHERE( (rmskkatax < 0.000001) .AND. (rmskkatax > -0.000001) ) 519 rmskkatax=1 520 rmskkatay=1 521 END WHERE 522 523 CALL lbc_lnk( rmskkatax(:,:), 'T', 1. ) ; CALL lbc_lnk( rmskkatay(:,:), 'T', 1. ) 524 525 IF(MAXVAL(rmskkatax) > 6.00001) THEN 526 WRITE(ctmp1,*) 'Problem in the data mask : maxval = ',MAXVAL(rmskkatax),' ( it is GT 6)' 527 CALL ctl_stop( ctmp1 ) 528 ENDIF 529 ENDIF 530 531 ! apply katabatic wind enhancement on grid T (before projection) 532 tauxt(:,:) = rmskkatax(:,:) * tauxt(:,:) 533 tauyt(:,:) = rmskkatay(:,:) * tauyt(:,:) 534 ! 535 ENDIF 536 ! 468 537 !Tau_x at U-point 469 538 DO ji = 1, jpim1 … … 486 555 ! --------------------------------- 487 556 ! 488 ! Sensible Heat : 489 qsb_oce(:,:) = rhoa * cp * Ch(:,:) * ( zsst(:,:) - flxnow(:,:,7) ) * dUnormt(:,:) 490 ! 491 ! Latent Heat : 492 evap(:,:) = rhoa * Ce(:,:) * ( qsatw(:,:) - flxnow(:,:,4) ) * dUnormt(:,:) 493 qla_oce(:,:) = Lv * evap(:,:) 494 ! 557 IF ( ln_2m ) THEN 558 !! 559 !! Values of temp. and hum. adjusted to 10m must be used instead of 2m values 560 !! Sensible Heat : ! right sign for ocean 561 qsb_oce(:,:) = rhoa*cp*Ch(:,:)*(zt_zu(:,:) - zsst(:,:))*dUnormt(:,:) 562 !! 563 !! Latent Heat : ! wrong sign for ocean 564 evap(:,:) = rhoa*Ce(:,:)*(qsatw(:,:) - zq_zu(:,:))*dUnormt(:,:) 565 !! 566 ELSE 567 !! 568 !! Sensible Heat : ! right sign for ocean 569 qsb_oce(:,:) = rhoa*cp*Ch(:,:)*(flxnow(:,:,7) - zsst(:,:))*dUnormt(:,:) 570 !! 571 !! Latent Heat : ! wrong sign for ocean 572 evap(:,:) = rhoa*Ce(:,:)*(qsatw(:,:) - flxnow(:,:,4))*dUnormt(:,:) 573 !! 574 END IF 575 !! 576 qla_oce(:,:) = Lv * evap(:,:) ! right sign for ocean 577 495 578 ! 496 579 ! II.3) Turbulent fluxes over ice … … 620 703 621 704 622 SUBROUTINE TURB_CORE ( zzu, T_0, T_a, q_sat, q_a, &623 & 705 SUBROUTINE TURB_CORE_1Z( zzu, T_0, T_a, q_sat, q_a, & 706 & dU , C_d, C_h , C_e ) 624 707 !!---------------------------------------------------------------------- 625 708 !! *** ROUTINE turb_core *** 626 709 !! 627 710 !! ** Purpose : Computes turbulent transfert coefficients of surface 628 !! fluxes according to Large & Yeager (2004)711 !! fluxes according to Large & Yeager (2004). 629 712 !! 630 713 !! ** Method : I N E R T I A L D I S S I P A T I O N M E T H O D … … 679 762 xlogt 680 763 681 INTEGER :: jkk,ji,jj,ii,ij,ilocu(2), jk ! doomy loop for iterations 682 INTEGER, PARAMETER :: nit = 3 ! number of iterations 683 REAL(wp), DIMENSION(jpi,jpj) :: dbg1,dbg2,dbg3,dbg4 684 REAL :: zXmax 764 INTEGER :: jk ! doomy loop for iterations 765 INTEGER, PARAMETER :: nit = 3 ! number of iterations 685 766 686 767 INTEGER, DIMENSION(jpi,jpj) :: & 687 & stab, &! stability test integer688 & stabit! stability within iterative loop689 690 REAL(wp), PARAMETER :: 691 & pi = 3.14159, &! Pi692 & grav = 9.8, &! gravity693 & kappa = 0.4! von Karman's constant768 stab, & ! stability test integer 769 stabit ! stability within iterative loop 770 771 REAL(wp), PARAMETER :: & 772 pi = 3.14159, & ! Pi 773 grav = 9.8, & ! gravity 774 kappa = 0.4 ! von Karman's constant 694 775 !!---------------------------------------------------------------------- 695 776 !! … … 807 888 C_e = Ce 808 889 !! 809 END SUBROUTINE TURB_CORE 890 END SUBROUTINE TURB_CORE_1Z 810 891 811 892 893 SUBROUTINE TURB_CORE_2Z(zt, zu, sst, T_zt, q_sat, q_zt, dU, Cd, Ch, Ce, T_zu, q_zu) 894 !!---------------------------------------------------------------------- 895 !! *** ROUTINE turb_core *** 896 !! 897 !! ** Purpose : Computes turbulent transfert coefficients of surface 898 !! fluxes according to Large & Yeager (2004). 899 !! 900 !! ** Method : I N E R T I A L D I S S I P A T I O N M E T H O D 901 !! Momentum, Latent and sensible heat exchange coefficients 902 !! Caution: this procedure should only be used in cases when air 903 !! temperature (T_air) and air specific humidity (q_air) are at 2m 904 !! whereas wind (dU) is at 10m. 905 !! 906 !! References : 907 !! Large & Yeager, 2004 : ??? 908 !! History : 909 !! ! XX-XX (??? ) Original code 910 !! 9.0 ! 05-08 (L. Brodeau) Optimisation 911 !!---------------------------------------------------------------------- 912 !! * Arguments 913 REAL(wp), INTENT( in ) :: & 914 zt, & ! height for T_zt and q_zt [m] 915 zu ! height for dU [m] 916 !! 917 REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: & 918 sst, & ! sea surface temperature [Kelvin] 919 T_zt, & ! potential air temperature [Kelvin] 920 q_sat, & ! sea surface specific humidity [kg/kg] 921 q_zt, & ! specific air humidity [kg/kg] 922 dU ! wind module |U(zu)-U(0)| [m/s] 923 !! 924 REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: & 925 Cd, Ch, Ce, & 926 T_zu, & ! air temp. shifted at zu [K] 927 q_zu ! spec. hum. shifted at zu [kg/kg] 928 929 !! * Local declarations 930 REAL(wp), DIMENSION(jpi,jpj) :: & 931 dU10, & ! dU [m/s] 932 dT, & ! air/sea temperature differeence [K] 933 dq, & ! air/sea humidity difference [K] 934 Cd_n10, & ! 10m neutral drag coefficient 935 Ce_n10, & ! 10m neutral latent coefficient 936 Ch_n10, & ! 10m neutral sensible coefficient 937 sqrt_Cd_n10, & ! root square of Cd_n10 938 sqrt_Cd, & ! root square of Cd 939 T_vpot_u, & ! virtual potential temperature [K] 940 T_star, & ! turbulent scale of tem. fluct. 941 q_star, & ! turbulent humidity of temp. fluct. 942 U_star, & ! turb. scale of velocity fluct. 943 L, & ! Monin-Obukov length [m] 944 zeta_u, & ! stability parameter at height zu 945 zeta_t, & ! stability parameter at height zt 946 U_n10, & ! neutral wind velocity at 10m [m] 947 xlogt, xct 948 !! 949 INTEGER, PARAMETER :: nb_itt = 3 ! number of itterations 950 !! 951 !! Some physical constants : 952 REAL(wp), PARAMETER :: & 953 grav = 9.8, & ! gravity 954 kappa = 0.4 ! von Karman's constant 955 !! 956 INTEGER :: j_itt 957 INTEGER, DIMENSION(jpi,jpj) :: & 958 stab ! 1st stability test integer 959 !!---------------------------------------------------------------------- 960 !! 961 !! ------------- 962 !! S T A R T 963 !! ------------- 964 !! 965 !! I.1 ) Preliminary stuffs 966 !! ------------------------ 967 !! 968 !! Initial air/sea differences 969 dU10 = max(0.5, dU) ; dT = T_zt - sst ; dq = q_zt - q_sat 970 !! 971 !! Neutral Drag Coefficient : 972 stab = 0.5 + sign(0.5,dT) ! stab = 1 if dT > 0 -> STABLE 973 Cd_n10 = 1E-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 974 sqrt_Cd_n10 = sqrt(Cd_n10) 975 Ce_n10 = 1E-3*( 34.6 * sqrt_Cd_n10 ) 976 Ch_n10 = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1 - stab)) 977 !! 978 !! I.2 ) Computing Neutral Drag Coefficient 979 !! ---------------------------------------- 980 !! Initializing transf. coeff. with their first guess neutral equivalents : 981 Cd = Cd_n10 ; Ce = Ce_n10 ; Ch = Ch_n10 ; sqrt_Cd = sqrt(Cd) 982 !! 983 !! Initializing z_u values with z_t values : 984 T_zu = T_zt ; q_zu = q_zt 985 !! 986 !! II ) Now starting iteration loop (IDM) 987 !! ---------------------------------------- 988 !! 989 DO j_itt=1, nb_itt 990 !! 991 !! Updating air/sea differences : 992 dT = T_zu - sst ; dq = q_zu - q_sat 993 !! 994 !! Updating virtual potential temperature at zu : 995 T_vpot_u = T_zu*(1. + 0.608*q_zu) 996 !! 997 !! Updating turbulent scales : (L & Y eq. (7)) 998 U_star = sqrt_Cd*dU10 ; T_star = Ch/sqrt_Cd*dT ; q_star = Ce/sqrt_Cd*dq 999 !! 1000 !! Estimate the Monin-Obukov length at height zu : 1001 L = (U_star*U_star) & 1002 & / (kappa*grav/T_vpot_u*(T_star*(1.+0.608*q_zu) + 0.608*T_zu*q_star)) 1003 !! 1004 !! Stability parameters : 1005 zeta_u = zu/L ; zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 1006 zeta_t = zt/L ; zeta_t = sign( min(abs(zeta_t),10.0), zeta_t ) 1007 !! 1008 !! Shifting the wind speed to 10m and neutral stability : (L & Y eq.(9a)) 1009 U_n10 = dU10/(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u))) 1010 !! 1011 !! Shifting temperature and humidity at zu : (L & Y eq. (9b-9c)) 1012 T_zu = T_zt - T_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t)) 1013 q_zu = q_zt - q_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t)) 1014 !! 1015 !! q_zu cannot have a negative value : forcing 0 1016 stab = 0.5 + sign(0.5,q_zu) ; q_zu = stab*q_zu 1017 !! 1018 !! Updating the neutral 10m transfer coefficients : 1019 Cd_n10 = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a) 1020 sqrt_Cd_n10 = sqrt(Cd_n10) 1021 Ce_n10 = 1E-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b) 1022 stab = 0.5 + sign(0.5,zeta_u) 1023 Ch_n10 = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! L & Y eq. (6c-6d) 1024 !! 1025 !! 1026 !! Shifting the neutral 10m transfer coefficients to (zu,zeta_u) : 1027 xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u)) 1028 Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd) 1029 !! 1030 xlogt = log(zu/10.) - psi_h(zeta_u) 1031 !! 1032 xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10 1033 Ch = Ch_n10*sqrt_Cd/sqrt_Cd_n10/xct 1034 !! 1035 xct = 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10 1036 Ce = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct 1037 !! 1038 !! 1039 END DO 1040 !! 1041 END SUBROUTINE TURB_CORE_2Z 1042 1043 FUNCTION psi_m(zta) !! Psis, L & Y eq. (8c), (8d), (8e) 1044 REAL(wp), PARAMETER :: pi = 3.14159 1045 REAL(wp), DIMENSION(jpi,jpj) :: psi_m 1046 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta 1047 REAL(wp), DIMENSION(jpi,jpj) :: X2, X, stabit 1048 X2 = sqrt(abs(1. - 16.*zta)) ; X2 = max(X2 , 1.0) ; X = sqrt(X2) 1049 stabit = 0.5 + sign(0.5,zta) 1050 psi_m = -5.*zta*stabit & ! Stable 1051 & + (1. - stabit)*(2*log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2) ! Unstable 1052 END FUNCTION psi_m 1053 1054 FUNCTION psi_h(zta) !! Psis, L & Y eq. (8c), (8d), (8e) 1055 REAL(wp), DIMENSION(jpi,jpj) :: psi_h 1056 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta 1057 REAL(wp), DIMENSION(jpi,jpj) :: X2, X, stabit 1058 X2 = sqrt(abs(1. - 16.*zta)) ; X2 = max(X2 , 1.) ; X = sqrt(X2) 1059 stabit = 0.5 + sign(0.5,zta) 1060 psi_h = -5.*zta*stabit & ! Stable 1061 & + (1. - stabit)*(2.*log( (1. + X2)/2. )) ! Unstable 1062 END FUNCTION psi_h 1063 1064 812 1065 SUBROUTINE flx_blk_albedo( palb , palcn , palbp , palcnp ) 813 1066 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
