Changeset 441 for LMDZ.3.3/branches/rel-LF/libf
- Timestamp:
- Jan 27, 2003, 11:26:36 AM (22 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90
r398 r441 57 57 & klon, iim, jjm, nisurf, knon, knindex, pctsrf, & 58 58 & rlon, rlat, cufi, cvfi,& 59 & debut, lafin, ok_veget, soil_model, nsoilmx, tsoil, &59 & debut, lafin, ok_veget, soil_model, nsoilmx, tsoil, qsol,& 60 60 & zlev, u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, & 61 61 & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, & 62 62 & precip_rain, precip_snow, sollw, sollwdown, swnet, swdown, & 63 63 & fder, taux, tauy, rugos, rugoro, & 64 & albedo, snow, qs ol, &64 & albedo, snow, qsurf, & 65 65 & tsurf, p1lay, ps, radsol, & 66 66 & ocean, npas, nexca, zmasq, & … … 171 171 character (len = 6) :: ocean 172 172 integer :: npas, nexca ! nombre et pas de temps couplage 173 real, dimension(klon), intent(INOUT) :: evap, snow, qs ol173 real, dimension(klon), intent(INOUT) :: evap, snow, qsurf 174 174 !! PB ajout pour soil 175 175 logical :: soil_model 176 176 integer :: nsoilmx 177 177 REAL, DIMENSION(klon, nsoilmx) :: tsoil 178 REAL, dimension(klon), intent(INOUT) :: qsol 178 179 REAL, dimension(klon) :: soilcap 179 180 REAL, dimension(klon) :: soilflux … … 193 194 integer, save :: error 194 195 integer :: ii, index 195 logical,save :: check = . false.196 logical,save :: check = .true. 196 197 real, dimension(klon):: cal, beta, dif_grnd, capsol 197 198 !!$PB real, parameter :: calice=1.0/(5.1444e+06*0.15), tau_gl=86400.*5. … … 200 201 real, dimension(klon):: alb_ice 201 202 real, dimension(klon):: tsurf_temp 202 real, dimension(klon):: qs ol_new203 real, dimension(klon):: qsurf_new 203 204 !! real, allocatable, dimension(:), save :: alb_neig_grid 204 205 real, dimension(klon):: alb_neig, alb_eau … … 309 310 & alb_new, z0_new) 310 311 ! 311 ! calcul snow et qs ol, hydrol adapté312 ! calcul snow et qsurf, hydrol adapté 312 313 ! 313 314 CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd) … … 323 324 CALL calcul_fluxs( klon, knon, nisurf, dtime, & 324 325 & tsurf, p1lay, cal, beta, tq_cdrag, ps, & 325 & precip_rain, precip_snow, snow, qs ol, &326 & precip_rain, precip_snow, snow, qsurf, & 326 327 & radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, & 327 328 & petAcoef, peqAcoef, petBcoef, peqBcoef, & … … 358 359 & evap, fluxsens, fluxlat, & 359 360 & tsol_rad, tsurf_new, alb_new, alblw, & 360 & emis_new, z0_new, dflux_l, dflux_s, qs ol_new)361 & emis_new, z0_new, dflux_l, dflux_s, qsurf_new) 361 362 362 363 ! … … 366 367 ! 367 368 ! mise a jour de l'humidite saturante calculee par ORCHIDEE 368 qs ol(1:knon) = qsol_new(1:knon)369 qsurf(1:knon) = qsurf_new(1:knon) 369 370 370 371 endif … … 416 417 call calcul_fluxs( klon, knon, nisurf, dtime, & 417 418 & tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, & 418 & precip_rain, precip_snow, snow, qs ol, &419 & precip_rain, precip_snow, snow, qsurf, & 419 420 & radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, & 420 421 & petAcoef, peqAcoef, petBcoef, peqBcoef, & … … 511 512 CALL calcul_fluxs( klon, knon, nisurf, dtime, & 512 513 & tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, & 513 & precip_rain, precip_snow, snow, qs ol, &514 & precip_rain, precip_snow, snow, qsurf, & 514 515 & radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, & 515 516 & petAcoef, peqAcoef, petBcoef, peqBcoef, & … … 581 582 call calcul_fluxs( klon, knon, nisurf, dtime, & 582 583 & tsurf, p1lay, cal, beta, tq_cdrag, ps, & 583 & precip_rain, precip_snow, snow, qs ol, &584 & precip_rain, precip_snow, snow, qsurf, & 584 585 & radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, & 585 586 & petAcoef, peqAcoef, petBcoef, peqBcoef, & … … 651 652 & knindex, rlon, rlat, cufi, cvfi, iim, jjm, pctsrf, & 652 653 & debut, lafin, ok_veget, & 653 & zlev, u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, &654 & plev, u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, & 654 655 & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, & 655 656 & precip_rain, precip_snow, lwdown, swnet, swdown, & … … 680 681 ! ok_veget logical: appel ou non au schema de surface continental 681 682 ! (si false calcul simplifie des fluxs sur les continents) 682 ! zlev hauteur de la premiere couche683 ! plev hauteur de la premiere couche (Pa) 683 684 ! u1_lay vitesse u 1ere couche 684 685 ! v1_lay vitesse v 1ere couche … … 715 716 ! emis_new emissivite 716 717 ! z0_new surface roughness 717 ! qsurf saturatedair moisture at surface718 ! qsurf air moisture at surface 718 719 719 720 ! Parametres d'entree … … 730 731 real, dimension(klon), intent(IN) :: rlon, rlat 731 732 real, dimension(klon), intent(IN) :: cufi, cvfi 732 real, dimension(klon), intent(IN) :: zlev733 real, dimension(klon), intent(IN) :: plev 733 734 real, dimension(klon), intent(IN) :: u1_lay, v1_lay 734 735 real, dimension(klon), intent(IN) :: temp_air, spechum … … 739 740 real, dimension(klon), intent(IN) :: precip_rain, precip_snow 740 741 real, dimension(klon), intent(IN) :: lwdown, swnet, swdown, ps 742 !IM cf. JP +++ 743 real, dimension(klon) :: swdown_vrai 744 !IM cf. JP --- 741 745 real, dimension(klon), intent(IN) :: tsurf, p1lay 742 746 real, dimension(klon), intent(IN) :: radsol … … 759 763 ! drapeaux controlant les appels dans SECHIBA 760 764 ! type(control_type), save :: control_in 765 ! Preserved albedo 766 !IM cf. JP +++ 767 real, allocatable, dimension(:), save :: albedo_keep, zlev 768 !IM cf. JP --- 761 769 ! coordonnees geographiques 762 770 real, allocatable, dimension(:,:), save :: lalo … … 794 802 795 803 #include "temps.inc" 804 #include "YOMCST.inc" 796 805 797 806 if (check) write(*,*)'Entree ', modname … … 803 812 if (debut) then 804 813 814 IF ( .NOT. allocated(albedo_keep)) THEN 815 ALLOCATE(albedo_keep(klon)) 816 ALLOCATE(zlev(klon)) 817 ENDIF 805 818 ! Pb de correspondances de grilles 806 819 allocate(ig(klon)) … … 990 1003 cdrag(1:knon) = tq_cdrag(1:knon) 991 1004 1005 !IM cf. JP +++ 1006 ! zlev(1:knon) = (100.*plev(1:knon))/((ps(1:knon)/287.05*temp_air(1:knon))*9.80665) 1007 zlev(1:knon) = (100.*plev(1:knon))/((ps(1:knon)/RD*temp_air(1:knon))*RG) 1008 !IM cf. JP --- 1009 992 1010 where(cdrag > 0.01) 993 1011 cdrag = 0.01 … … 1008 1026 & tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, & 1009 1027 & lon_scat, lat_scat) 1028 1029 !IM cf. JP +++ 1030 albedo_keep(:) = (albedo_out(:,1)+albedo_out(:,2))/2. 1031 !IM cf. JP --- 1032 1010 1033 endif 1034 1035 !IM cf. JP +++ 1036 swdown_vrai(:) = swnet(:)/(1. - albedo_keep(:)) 1037 !IM cf. JP --- 1011 1038 1012 1039 call intersurf_main (itime+itau_phy, iim, jjm+1, knon, ktindex, dtime, & … … 1015 1042 & zlev, u1_lay, v1_lay, spechum, temp_air, epot_air, ccanopy, & 1016 1043 & cdrag, petA_orc, peqA_orc, petB_orc, peqB_orc, & 1017 & precip_rain, precip_snow, lwdown, swnet, swdown, ps, & 1044 !IM cf. JP +++ 1045 & precip_rain, precip_snow, lwdown, swnet, swdown_vrai, ps, & 1046 !IM cf. JP --- 1018 1047 & evap, fluxsens, fluxlat, coastalflow, riverflow, & 1019 1048 & tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, & 1020 1049 & lon_scat, lat_scat) 1050 1051 !IM cf. JP +++ 1052 albedo_keep(:) = (albedo_out(:,1)+albedo_out(:,2))/2. 1053 !IM cf. JP --- 1021 1054 1022 1055 bidule=0. … … 1927 1960 character (len = 20),save :: fich ='limit.nc' 1928 1961 logical,save :: newlmt = .false. 1929 logical,save :: check = . FALSE.1962 logical,save :: check = .true. 1930 1963 ! Champs lus dans le fichier de CL 1931 1964 real, allocatable , save, dimension(:) :: alb_lu, rug_lu … … 1951 1984 if (check) write(*,*)modname,':: jour_lu_sur, deja_lu_sur', jour_lu_sur, deja_lu_sur 1952 1985 if (check) write(*,*)modname,':: itime, lmt_pas', itime, lmt_pas 1953 call flush(6)1986 if (check) call flush(6) 1954 1987 1955 1988 ! Tester d'abord si c'est le moment de lire le fichier … … 1959 1992 ! 1960 1993 fich = trim(fich) 1994 IF (check) WRITE(*,*)modname,' ouverture fichier ',fich 1995 if (check) CALL flush(6) 1961 1996 ierr = NF_OPEN (fich, NF_NOWRITE,nid) 1962 1997 if (ierr.NE.NF_NOERR) then … … 2033 2068 SUBROUTINE calcul_fluxs( klon, knon, nisurf, dtime, & 2034 2069 & tsurf, p1lay, cal, beta, coef1lay, ps, & 2035 & precip_rain, precip_snow, snow, qs ol, &2070 & precip_rain, precip_snow, snow, qsurf, & 2036 2071 & radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, & 2037 2072 & petAcoef, peqAcoef, petBcoef, peqBcoef, & … … 2055 2090 ! precip_snow precipitations solides 2056 2091 ! snow champs hauteur de neige 2057 ! qsol humidite du sol2058 2092 ! runoff runoff en cas de trop plein 2059 2093 ! petAcoef coeff. A de la resolution de la CL pour t … … 2066 2100 ! output: 2067 2101 ! tsurf_new temperature au sol 2102 ! qsurf humidite de l'air au dessus du sol 2068 2103 ! fluxsens flux de chaleur sensible 2069 2104 ! fluxlat flux de chaleur latente … … 2086 2121 real, dimension(klon), intent(IN) :: radsol, dif_grnd 2087 2122 real, dimension(klon), intent(IN) :: t1lay, u1lay, v1lay 2088 real, dimension(klon), intent(INOUT) :: snow, qs ol2123 real, dimension(klon), intent(INOUT) :: snow, qsurf 2089 2124 2090 2125 ! Parametres sorties … … 2102 2137 real :: bilan_f, fq_fonte 2103 2138 REAL :: subli, fsno 2139 REAL :: qsat_new, q1_new 2104 2140 real, parameter :: t_grnd = 271.35, t_coup = 273.15 2105 2141 !! PB temporaire en attendant mieux pour le modele de neige 2106 2142 REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15) 2107 2143 ! 2108 logical, save :: check = . FALSE.2144 logical, save :: check = .true. 2109 2145 character (len = 20) :: modname = 'calcul_fluxs' 2110 2146 logical, save :: fonte_neige = .false. … … 2114 2150 2115 2151 if (check) write(*,*)'Entree ', modname,' surface = ',nisurf 2152 2153 IF (check) THEN 2154 WRITE(*,*)' radsol (min, max)' & 2155 & , MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon)) 2156 CALL flush(6) 2157 ENDIF 2116 2158 2117 2159 if (size(coastalflow) /= knon .AND. nisurf == is_ter) then … … 2135 2177 !!$ where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime 2136 2178 !!$ endif 2137 IF (nisurf /= is_ter) qsol = max_eau_sol2179 !!$ IF (nisurf /= is_ter) qsol = max_eau_sol 2138 2180 2139 2181 … … 2233 2275 dflux_s(i) = zx_nh(i) 2234 2276 dflux_l(i) = (zx_sl(i) * zx_nq(i)) 2235 2277 ! Nouvelle valeure de l'humidite au dessus du sol 2278 qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i) 2279 q1_new = peqAcoef(i) - peqBcoef(i)*evap(i)*dtime 2280 qsurf(i)=q1_new*(1.-beta(i)) + beta(i)*qsat_new 2236 2281 ! 2237 2282 ! en cas de fonte de neige … … 2412 2457 ! precip_snow precipitations solides 2413 2458 ! snow champs hauteur de neige 2414 ! qsol h umidite dusol2459 ! qsol hauteur d'eau contenu dans le sol 2415 2460 ! runoff runoff en cas de trop plein 2416 2461 ! petAcoef coeff. A de la resolution de la CL pour t … … 2459 2504 real :: bilan_f, fq_fonte 2460 2505 REAL :: subli, fsno 2506 real, dimension(klon) :: bil_eau_s 2461 2507 real, parameter :: t_grnd = 271.35, t_coup = 273.15 2462 2508 !! PB temporaire en attendant mieux pour le modele de neige … … 2473 2519 2474 2520 ! Initialisations 2521 bil_eau_s(:) = 0. 2475 2522 DO i = 1, knon 2476 2523 zx_pkh(i) = (ps(i)/ps(i))**RKAPPA … … 2534 2581 WHERE (precip_snow > 0.) snow = snow + (precip_snow * dtime) 2535 2582 WHERE (evap > 0 ) snow = MAX(0.0, snow - (evap * dtime)) 2536 qsol = qsol+ (precip_rain - evap) * dtime2583 bil_eau_s = bil_eau_s + (precip_rain - evap) * dtime 2537 2584 ! 2538 2585 ! Y'a-t-il fonte de neige? … … 2544 2591 fq_fonte = MIN( MAX((tsurf_new(i)-RTT )/chasno,0.0),snow(i)) 2545 2592 snow(i) = max(0., snow(i) - fq_fonte) 2546 qsol(i) = qsol(i) + fq_fonte2593 bil_eau_s(i) = bil_eau_s(i) + fq_fonte 2547 2594 tsurf_new(i) = tsurf_new(i) - fq_fonte * chasno 2548 2595 IF (nisurf == is_sic .OR. nisurf == is_lic ) tsurf_new(i) = RTT -1.8 … … 2564 2611 !!$ fq_fonte = bilan_f / zx_sl(i) 2565 2612 endif 2566 IF (nisurf == is_ter) & 2567 & run_off(i) = run_off(i) + MAX(qsol(i) - max_eau_sol, 0.0) 2568 qsol(i) = MIN(qsol(i), max_eau_sol) 2613 2614 IF (nisurf == is_ter) then 2615 qsol(i) = qsol(i) + bil_eau_s(i) 2616 run_off(i) = run_off(i) + MAX(qsol(i) - max_eau_sol, 0.0) 2617 qsol(i) = MIN(qsol(i), max_eau_sol) 2618 else 2619 run_off(i) = run_off(i) + MAX(bil_eau_s(i), 0.0) 2620 endif 2569 2621 enddo 2570 2622
Note: See TracChangeset
for help on using the changeset viewer.