--- trunk/libf/phylmd/interface_surf.f90 2011/10/07 13:11:58 53 +++ trunk/libf/phylmd/Interface_surf/interface_surf.f90 2011/12/06 15:07:04 54 @@ -1,21 +1,10 @@ MODULE interface_surf ! From phylmd/interface_surf.F90, version 1.8 2005/05/25 13:10:09 - - ! Ce module regroupe toutes les routines gérant l'interface entre le - ! modèle atmosphérique et les modèles de surface (sols continentaux, - ! océans, glaces). Les routines sont les suivantes. "interfsurf_hq" : - ! routine d'aiguillage vers les interfaces avec les différents - ! modèles de surface ; "interfoce_*" : routines d'interface proprement - ! dites. - ! L. Fairhead, LMD, february 2000 IMPLICIT none - PRIVATE - PUBLIC :: interfsurf_hq - ! run_off ruissellement total REAL, ALLOCATABLE, DIMENSION(:), SAVE :: run_off, run_off_lic real, allocatable, dimension(:), save :: coastalflow, riverflow @@ -24,1638 +13,4 @@ ! pour simuler la fonte des glaciers antarctiques REAL, save :: tau_calv -CONTAINS - - SUBROUTINE interfsurf_hq(itime, dtime, date0, jour, rmu0, klon, iim, jjm, & - nisurf, knon, knindex, pctsrf, rlon, rlat, cufi, cvfi, debut, lafin, & - ok_veget, soil_model, nsoilmx, tsoil, qsol, zlev, u1_lay, v1_lay, & - temp_air, spechum, epot_air, ccanopy, tq_cdrag, petAcoef, peqAcoef, & - petBcoef, peqBcoef, precip_rain, precip_snow, sollw, sollwdown, swnet, & - swdown, fder, taux, tauy, windsp, rugos, rugoro, albedo, snow, qsurf, & - tsurf, p1lay, ps, radsol, ocean, npas, nexca, zmasq, evap, fluxsens, & - fluxlat, dflux_l, dflux_s, tsol_rad, tsurf_new, alb_new, alblw, & - emis_new, z0_new, pctsrf_new, agesno, fqcalving, ffonte, & - run_off_lic_0, flux_o, flux_g, tslab, seaice) - - ! Cette routine sert d'aiguillage entre l'atmosphère et la surface - ! en général (sols continentaux, océans, glaces) pour les flux de - ! chaleur et d'humidité. - ! En pratique l'interface se fait entre la couche limite du modèle - ! atmosphérique ("clmain.F") et les routines de surface - ! ("sechiba", "oasis"...). - - ! L.Fairhead 02/2000 - - use abort_gcm_m, only: abort_gcm - use gath_cpl, only: gath2cpl - use indicesol - use SUPHEC_M - use albsno_m, only: albsno - - ! Parametres d'entree - ! input: - ! klon nombre total de points de grille - ! iim, jjm nbres de pts de grille - ! dtime pas de temps de la physique (en s) - ! date0 jour initial - ! jour jour dans l'annee en cours, - ! rmu0 cosinus de l'angle solaire zenithal - ! nexca pas de temps couplage - ! nisurf index de la surface a traiter (1 = sol continental) - ! knon nombre de points de la surface a traiter - ! knindex index des points de la surface a traiter - ! pctsrf tableau des pourcentages de surface de chaque maille - ! rlon longitudes - ! rlat latitudes - ! cufi, cvfi resolution des mailles en x et y (m) - ! debut logical: 1er appel a la physique - ! lafin logical: dernier appel a la physique - ! ok_veget logical: appel ou non au schema de surface continental - ! (si false calcul simplifie des fluxs sur les continents) - ! zlev hauteur de la premiere couche - ! u1_lay vitesse u 1ere couche - ! v1_lay vitesse v 1ere couche - ! temp_air temperature de l'air 1ere couche - ! spechum humidite specifique 1ere couche - ! epot_air temp potentielle de l'air - ! ccanopy concentration CO2 canopee - ! tq_cdrag cdrag - ! petAcoef coeff. A de la resolution de la CL pour t - ! peqAcoef coeff. A de la resolution de la CL pour q - ! petBcoef coeff. B de la resolution de la CL pour t - ! peqBcoef coeff. B de la resolution de la CL pour q - ! precip_rain precipitation liquide - ! precip_snow precipitation solide - ! sollw flux IR net a la surface - ! sollwdown flux IR descendant a la surface - ! swnet flux solaire net - ! swdown flux solaire entrant a la surface - ! albedo albedo de la surface - ! tsurf temperature de surface - ! tslab temperature slab ocean - ! pctsrf_slab pourcentages (0-1) des sous-surfaces dans le slab - ! tmp_pctsrf_slab = pctsrf_slab - ! p1lay pression 1er niveau (milieu de couche) - ! ps pression au sol - ! radsol rayonnement net aus sol (LW + SW) - ! ocean type d'ocean utilise ("force" ou "slab" mais pas "couple") - ! fder derivee des flux (pour le couplage) - ! taux, tauy tension de vents - ! windsp module du vent a 10m - ! rugos rugosite - ! zmasq masque terre/ocean - ! rugoro rugosite orographique - ! run_off_lic_0 runoff glacier du pas de temps precedent - integer, intent(IN) :: itime ! numero du pas de temps - integer, intent(IN) :: iim, jjm - integer, intent(IN) :: klon - real, intent(IN) :: dtime - real, intent(IN) :: date0 - integer, intent(IN) :: jour - real, intent(IN) :: rmu0(klon) - integer, intent(IN) :: nisurf - integer, intent(IN) :: knon - integer, dimension(klon), intent(in) :: knindex - real, dimension(klon, nbsrf), intent(IN) :: pctsrf - logical, intent(IN) :: debut, lafin, ok_veget - real, dimension(klon), intent(IN) :: rlon, rlat - real, dimension(klon), intent(IN) :: cufi, cvfi - real, dimension(klon), intent(INOUT) :: tq_cdrag - real, dimension(klon), intent(IN) :: zlev - real, dimension(klon), intent(IN) :: u1_lay, v1_lay - real, dimension(klon), intent(IN) :: temp_air, spechum - real, dimension(klon), intent(IN) :: epot_air, ccanopy - real, dimension(klon), intent(IN) :: petAcoef, peqAcoef - real, dimension(klon), intent(IN) :: petBcoef, peqBcoef - real, dimension(klon), intent(IN) :: precip_rain, precip_snow - real, dimension(klon), intent(IN) :: sollw, sollwdown, swnet, swdown - real, dimension(klon), intent(IN) :: ps, albedo - real, dimension(klon), intent(IN) :: tsurf, p1lay - !IM: "slab" ocean - real, dimension(klon), intent(INOUT) :: tslab - real, allocatable, dimension(:), save :: tmp_tslab - real, dimension(klon), intent(OUT) :: flux_o, flux_g - real, dimension(klon), intent(INOUT) :: seaice ! glace de mer (kg/m2) - REAL, DIMENSION(klon), INTENT(INOUT) :: radsol, fder - real, dimension(klon), intent(IN) :: zmasq - real, dimension(klon), intent(IN) :: taux, tauy, rugos, rugoro - real, dimension(klon), intent(IN) :: windsp - character(len=*), intent(IN):: ocean - integer :: npas, nexca ! nombre et pas de temps couplage - real, dimension(klon), intent(INOUT) :: evap, snow, qsurf - !! PB ajout pour soil - logical, intent(in):: soil_model - integer :: nsoilmx - REAL, DIMENSION(klon, nsoilmx) :: tsoil - REAL, dimension(klon), intent(INOUT) :: qsol - REAL, dimension(klon) :: soilcap - REAL, dimension(klon) :: soilflux - - ! Parametres de sortie - ! output: - ! evap evaporation totale - ! fluxsens flux de chaleur sensible - ! fluxlat flux de chaleur latente - ! tsol_rad - ! tsurf_new temperature au sol - ! alb_new albedo - ! emis_new emissivite - ! z0_new surface roughness - ! pctsrf_new nouvelle repartition des surfaces - real, dimension(klon), intent(OUT):: fluxsens, fluxlat - real, dimension(klon), intent(OUT):: tsol_rad, tsurf_new, alb_new - real, dimension(klon), intent(OUT):: alblw - real, dimension(klon), intent(OUT):: emis_new, z0_new - real, dimension(klon), intent(OUT):: dflux_l, dflux_s - real, dimension(klon, nbsrf), intent(OUT) :: pctsrf_new - real, dimension(klon), intent(INOUT):: agesno - real, dimension(klon), intent(INOUT):: run_off_lic_0 - - ! Flux thermique utiliser pour fondre la neige - !jld a rajouter real, dimension(klon), intent(INOUT):: ffonte - real, dimension(klon), intent(INOUT):: ffonte - ! Flux d'eau "perdue" par la surface et nécessaire pour que limiter la - ! hauteur de neige, en kg/m2/s - !jld a rajouter real, dimension(klon), intent(INOUT):: fqcalving - real, dimension(klon), intent(INOUT):: fqcalving - !IM: "slab" ocean - Local - real, parameter :: t_grnd=271.35 - real, dimension(klon) :: zx_sl - integer i - real, allocatable, dimension(:), save :: tmp_flux_o, tmp_flux_g - real, allocatable, dimension(:), save :: tmp_radsol - real, allocatable, dimension(:, :), save :: tmp_pctsrf_slab - real, allocatable, dimension(:), save :: tmp_seaice - - ! Local - character (len = 20), save :: modname = 'interfsurf_hq' - character (len = 80) :: abort_message - logical, save :: first_call = .true. - integer, save :: error - integer :: ii - logical, save :: check = .false. - real, dimension(klon):: cal, beta, dif_grnd, capsol - real, parameter :: calice=1.0/(5.1444e+06*0.15), tau_gl=86400.*5. - real, parameter :: calsno=1./(2.3867e+06*.15) - real, dimension(klon):: tsurf_temp - real, dimension(klon):: alb_neig, alb_eau - real, DIMENSION(klon):: zfra - logical :: cumul = .false. - INTEGER, dimension(1) :: iloc - real, dimension(klon):: fder_prev - REAL, dimension(klon) :: bidule - - !------------------------------------------------------------- - - if (check) write(*, *) 'Entree ', modname - - ! On doit commencer par appeler les schemas de surfaces continentales - ! car l'ocean a besoin du ruissellement qui est y calcule - - if (first_call) then - call conf_interface(tau_calv) - if (nisurf /= is_ter .and. klon > 1) then - write(*, *)' *** Warning ***' - write(*, *)' nisurf = ', nisurf, ' /= is_ter = ', is_ter - write(*, *)'or on doit commencer par les surfaces continentales' - abort_message='voir ci-dessus' - call abort_gcm(modname, abort_message, 1) - endif - if (ocean /= 'slab' .and. ocean /= 'force') then - write(*, *)' *** Warning ***' - write(*, *)'Option couplage pour l''ocean = ', ocean - abort_message='option pour l''ocean non valable' - call abort_gcm(modname, abort_message, 1) - endif - if ( is_oce > is_sic ) then - write(*, *)' *** Warning ***' - write(*, *)' Pour des raisons de sequencement dans le code' - write(*, *)' l''ocean doit etre traite avant la banquise' - write(*, *)' or is_oce = ', is_oce, '> is_sic = ', is_sic - abort_message='voir ci-dessus' - call abort_gcm(modname, abort_message, 1) - endif - endif - first_call = .false. - - ! Initialisations diverses - - ffonte(1:knon)=0. - fqcalving(1:knon)=0. - - cal = 999999. ; beta = 999999. ; dif_grnd = 999999. ; capsol = 999999. - alb_new = 999999. ; z0_new = 999999. ; alb_neig = 999999. - tsurf_new = 999999. - alblw = 999999. - - !IM: "slab" ocean; initialisations - flux_o = 0. - flux_g = 0. - - if (.not. allocated(tmp_flux_o)) then - allocate(tmp_flux_o(klon), stat = error) - DO i=1, knon - tmp_flux_o(knindex(i))=flux_o(i) - ENDDO - if (error /= 0) then - abort_message='Pb allocation tmp_flux_o' - call abort_gcm(modname, abort_message, 1) - endif - endif - if (.not. allocated(tmp_flux_g)) then - allocate(tmp_flux_g(klon), stat = error) - DO i=1, knon - tmp_flux_g(knindex(i))=flux_g(i) - ENDDO - if (error /= 0) then - abort_message='Pb allocation tmp_flux_g' - call abort_gcm(modname, abort_message, 1) - endif - endif - if (.not. allocated(tmp_radsol)) then - allocate(tmp_radsol(klon), stat = error) - if (error /= 0) then - abort_message='Pb allocation tmp_radsol' - call abort_gcm(modname, abort_message, 1) - endif - endif - DO i=1, knon - tmp_radsol(knindex(i))=radsol(i) - ENDDO - if (.not. allocated(tmp_pctsrf_slab)) then - allocate(tmp_pctsrf_slab(klon, nbsrf), stat = error) - if (error /= 0) then - abort_message='Pb allocation tmp_pctsrf_slab' - call abort_gcm(modname, abort_message, 1) - endif - DO i=1, klon - tmp_pctsrf_slab(i, 1:nbsrf)=pctsrf(i, 1:nbsrf) - ENDDO - endif - - if (.not. allocated(tmp_seaice)) then - allocate(tmp_seaice(klon), stat = error) - if (error /= 0) then - abort_message='Pb allocation tmp_seaice' - call abort_gcm(modname, abort_message, 1) - endif - DO i=1, klon - tmp_seaice(i)=seaice(i) - ENDDO - endif - - if (.not. allocated(tmp_tslab)) then - allocate(tmp_tslab(klon), stat = error) - if (error /= 0) then - abort_message='Pb allocation tmp_tslab' - call abort_gcm(modname, abort_message, 1) - endif - endif - DO i=1, klon - tmp_tslab(i)=tslab(i) - ENDDO - - ! Aiguillage vers les differents schemas de surface - - if (nisurf == is_ter) then - - ! Surface "terre" appel a l'interface avec les sols continentaux - - ! allocation du run-off - if (.not. allocated(coastalflow)) then - allocate(coastalflow(knon), stat = error) - if (error /= 0) then - abort_message='Pb allocation coastalflow' - call abort_gcm(modname, abort_message, 1) - endif - allocate(riverflow(knon), stat = error) - if (error /= 0) then - abort_message='Pb allocation riverflow' - call abort_gcm(modname, abort_message, 1) - endif - allocate(run_off(knon), stat = error) - if (error /= 0) then - abort_message='Pb allocation run_off' - call abort_gcm(modname, abort_message, 1) - endif - !cym - run_off=0.0 - !cym - -!!$PB - ALLOCATE (tmp_rriv(iim, jjm+1), stat=error) - if (error /= 0) then - abort_message='Pb allocation tmp_rriv' - call abort_gcm(modname, abort_message, 1) - endif - ALLOCATE (tmp_rcoa(iim, jjm+1), stat=error) - if (error /= 0) then - abort_message='Pb allocation tmp_rcoa' - call abort_gcm(modname, abort_message, 1) - endif - ALLOCATE (tmp_rlic(iim, jjm+1), stat=error) - if (error /= 0) then - abort_message='Pb allocation tmp_rlic' - call abort_gcm(modname, abort_message, 1) - endif - tmp_rriv = 0.0 - tmp_rcoa = 0.0 - tmp_rlic = 0.0 - -!!$ - else if (size(coastalflow) /= knon) then - write(*, *)'Bizarre, le nombre de points continentaux' - write(*, *)'a change entre deux appels. J''arrete ...' - abort_message='voir ci-dessus' - call abort_gcm(modname, abort_message, 1) - endif - coastalflow = 0. - riverflow = 0. - - ! Calcul age de la neige - - if (.not. ok_veget) then - ! calcul albedo: lecture albedo fichier boundary conditions - ! puis ajout albedo neige - call interfsur_lim(itime, dtime, jour, klon, nisurf, knon, knindex, & - debut, alb_new, z0_new) - - ! calcul snow et qsurf, hydrol adapté - CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd) - - IF (soil_model) THEN - CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil, soilcap, & - soilflux) - cal(1:knon) = RCPD / soilcap(1:knon) - radsol(1:knon) = radsol(1:knon) + soilflux(1:knon) - ELSE - cal = RCPD * capsol - ENDIF - CALL calcul_fluxs( klon, knon, nisurf, dtime, & - tsurf, p1lay, cal, beta, tq_cdrag, ps, & - precip_rain, precip_snow, snow, qsurf, & - radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, & - petAcoef, peqAcoef, petBcoef, peqBcoef, & - tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) - - CALL fonte_neige( klon, knon, nisurf, dtime, & - tsurf, p1lay, cal, beta, tq_cdrag, ps, & - precip_rain, precip_snow, snow, qsol, & - radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, & - petAcoef, peqAcoef, petBcoef, peqBcoef, & - tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & - fqcalving, ffonte, run_off_lic_0) - - call albsno(klon, knon, dtime, agesno, alb_neig, precip_snow) - where (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0. - zfra(1:knon) = max(0.0, min(1.0, snow(1:knon)/(snow(1:knon)+10.0))) - alb_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + & - alb_new(1 : knon)*(1.0-zfra(1:knon)) - z0_new = sqrt(z0_new**2+rugoro**2) - alblw(1 : knon) = alb_new(1 : knon) - endif - - ! Remplissage des pourcentages de surface - pctsrf_new(:, nisurf) = pctsrf(:, nisurf) - else if (nisurf == is_oce) then - ! Surface "ocean" appel a l'interface avec l'ocean - if (ocean == 'slab') then - tsurf_new = tsurf - pctsrf_new = tmp_pctsrf_slab - else - ! lecture conditions limites - call interfoce_lim(itime, dtime, jour, klon, nisurf, knon, knindex, & - debut, tsurf_new, pctsrf_new) - endif - - tsurf_temp = tsurf_new - cal = 0. - beta = 1. - dif_grnd = 0. - alb_neig = 0. - agesno = 0. - - call calcul_fluxs( klon, knon, nisurf, dtime, & - tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, & - precip_rain, precip_snow, snow, qsurf, & - radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, & - petAcoef, peqAcoef, petBcoef, peqBcoef, & - tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) - - fder_prev = fder - fder = fder_prev + dflux_s + dflux_l - - iloc = maxloc(fder(1:klon)) - if (check.and.fder(iloc(1))> 0.) then - WRITE(*, *)'**** Debug fder****' - WRITE(*, *)'max fder(', iloc(1), ') = ', fder(iloc(1)) - WRITE(*, *)'fder_prev, dflux_s, dflux_l', fder_prev(iloc(1)), & - dflux_s(iloc(1)), dflux_l(iloc(1)) - endif - - !IM: flux ocean-atmosphere utile pour le "slab" ocean - DO i=1, knon - zx_sl(i) = RLVTT - if (tsurf_new(i) .LT. RTT) zx_sl(i) = RLSTT - flux_o(i) = fluxsens(i)-evap(i)*zx_sl(i) - tmp_flux_o(knindex(i)) = flux_o(i) - tmp_radsol(knindex(i))=radsol(i) - ENDDO - - ! 2eme appel a interfoce pour le cumul des champs (en particulier - ! fluxsens et fluxlat calcules dans calcul_fluxs) - - if (ocean == 'slab ') then - seaice=tmp_seaice - cumul = .true. - call interfoce_slab(klon, debut, itime, dtime, jour, & - tmp_radsol, tmp_flux_o, tmp_flux_g, pctsrf, & - tslab, seaice, pctsrf_new) - - tmp_pctsrf_slab=pctsrf_new - DO i=1, knon - tsurf_new(i)=tslab(knindex(i)) - ENDDO - endif - - ! calcul albedo - if ( minval(rmu0) == maxval(rmu0) .and. minval(rmu0) == -999.999 ) then - CALL alboc(FLOAT(jour), rlat, alb_eau) - else ! cycle diurne - CALL alboc_cd(rmu0, alb_eau) - endif - DO ii =1, knon - alb_new(ii) = alb_eau(knindex(ii)) - enddo - - z0_new = sqrt(rugos**2 + rugoro**2) - alblw(1:knon) = alb_new(1:knon) - else if (nisurf == is_sic) then - if (check) write(*, *)'sea ice, nisurf = ', nisurf - - ! Surface "glace de mer" appel a l'interface avec l'ocean - - - if (ocean == 'slab ') then - pctsrf_new=tmp_pctsrf_slab - - DO ii = 1, knon - tsurf_new(ii) = tsurf(ii) - IF (pctsrf_new(knindex(ii), nisurf) < EPSFRA) then - snow(ii) = 0.0 - tsurf_new(ii) = RTT - 1.8 - IF (soil_model) tsoil(ii, :) = RTT -1.8 - ENDIF - ENDDO - - CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd) - - IF (soil_model) THEN - CALL soil(dtime, nisurf, knon, snow, tsurf_new, tsoil, soilcap, soilflux) - cal(1:knon) = RCPD / soilcap(1:knon) - radsol(1:knon) = radsol(1:knon) + soilflux(1:knon) - ELSE - dif_grnd = 1.0 / tau_gl - cal = RCPD * calice - WHERE (snow > 0.0) cal = RCPD * calsno - ENDIF - tsurf_temp = tsurf_new - beta = 1.0 - - ELSE - ! ! lecture conditions limites - CALL interfoce_lim(itime, dtime, jour, & - klon, nisurf, knon, knindex, & - debut, & - tsurf_new, pctsrf_new) - - !IM cf LF - DO ii = 1, knon - tsurf_new(ii) = tsurf(ii) - !IMbad IF (pctsrf_new(ii, nisurf) < EPSFRA) then - IF (pctsrf_new(knindex(ii), nisurf) < EPSFRA) then - snow(ii) = 0.0 - !IM cf LF/JLD tsurf(ii) = RTT - 1.8 - tsurf_new(ii) = RTT - 1.8 - IF (soil_model) tsoil(ii, :) = RTT -1.8 - endif - enddo - - CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd) - - IF (soil_model) THEN - !IM cf LF/JLD CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil, soilcap, soilflux) - CALL soil(dtime, nisurf, knon, snow, tsurf_new, tsoil, soilcap, soilflux) - cal(1:knon) = RCPD / soilcap(1:knon) - radsol(1:knon) = radsol(1:knon) + soilflux(1:knon) - dif_grnd = 0. - ELSE - dif_grnd = 1.0 / tau_gl - cal = RCPD * calice - WHERE (snow > 0.0) cal = RCPD * calsno - ENDIF - !IMbadtsurf_temp = tsurf - tsurf_temp = tsurf_new - beta = 1.0 - ENDIF - - CALL calcul_fluxs( klon, knon, nisurf, dtime, & - tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, & - precip_rain, precip_snow, snow, qsurf, & - radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, & - petAcoef, peqAcoef, petBcoef, peqBcoef, & - tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) - - !IM: flux entre l'ocean et la glace de mer pour le "slab" ocean - DO i = 1, knon - flux_g(i) = 0.0 - - !IM: faire dependre le coefficient de conduction de la glace de mer - ! de l'epaisseur de la glace de mer, dans l'hypothese ou le coeff. - ! actuel correspond a 3m de glace de mer, cf. L.Li - - ! IF(1.EQ.0) THEN - ! IF(siceh(i).GT.0.) THEN - ! new_dif_grnd(i) = dif_grnd(i)*3./siceh(i) - ! ELSE - ! new_dif_grnd(i) = 0. - ! ENDIF - ! ENDIF !(1.EQ.0) THEN - - IF (cal(i).GT.1.0e-15) flux_g(i)=(tsurf_new(i)-t_grnd) & - * dif_grnd(i) *RCPD/cal(i) - ! & * new_dif_grnd(i) *RCPD/cal(i) - tmp_flux_g(knindex(i))=flux_g(i) - tmp_radsol(knindex(i))=radsol(i) - ENDDO - - CALL fonte_neige( klon, knon, nisurf, dtime, & - tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, & - precip_rain, precip_snow, snow, qsol, & - radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, & - petAcoef, peqAcoef, petBcoef, peqBcoef, & - tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & - fqcalving, ffonte, run_off_lic_0) - - ! calcul albedo - - CALL albsno(klon, knon, dtime, agesno, alb_neig, precip_snow) - WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0. - zfra(1:knon) = MAX(0.0, MIN(1.0, snow(1:knon)/(snow(1:knon)+10.0))) - alb_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + & - 0.6 * (1.0-zfra(1:knon)) - - fder_prev = fder - fder = fder_prev + dflux_s + dflux_l - - iloc = maxloc(fder(1:klon)) - if (check.and.fder(iloc(1))> 0.) then - WRITE(*, *)'**** Debug fder ****' - WRITE(*, *)'max fder(', iloc(1), ') = ', fder(iloc(1)) - WRITE(*, *)'fder_prev, dflux_s, dflux_l', fder_prev(iloc(1)), & - dflux_s(iloc(1)), dflux_l(iloc(1)) - endif - - - ! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean - - z0_new = 0.002 - z0_new = SQRT(z0_new**2+rugoro**2) - alblw(1:knon) = alb_new(1:knon) - - else if (nisurf == is_lic) then - - if (check) write(*, *)'glacier, nisurf = ', nisurf - - if (.not. allocated(run_off_lic)) then - allocate(run_off_lic(knon), stat = error) - if (error /= 0) then - abort_message='Pb allocation run_off_lic' - call abort_gcm(modname, abort_message, 1) - endif - run_off_lic = 0. - endif - - ! Surface "glacier continentaux" appel a l'interface avec le sol - - IF (soil_model) THEN - CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil, soilcap, soilflux) - cal(1:knon) = RCPD / soilcap(1:knon) - radsol(1:knon) = radsol(1:knon) + soilflux(1:knon) - ELSE - cal = RCPD * calice - WHERE (snow > 0.0) cal = RCPD * calsno - ENDIF - beta = 1.0 - dif_grnd = 0.0 - - call calcul_fluxs( klon, knon, nisurf, dtime, & - tsurf, p1lay, cal, beta, tq_cdrag, ps, & - precip_rain, precip_snow, snow, qsurf, & - radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, & - petAcoef, peqAcoef, petBcoef, peqBcoef, & - tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) - - call fonte_neige( klon, knon, nisurf, dtime, & - tsurf, p1lay, cal, beta, tq_cdrag, ps, & - precip_rain, precip_snow, snow, qsol, & - radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, & - petAcoef, peqAcoef, petBcoef, peqBcoef, & - tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & - fqcalving, ffonte, run_off_lic_0) - - ! passage du run-off des glaciers calcule dans fonte_neige au coupleur - bidule=0. - bidule(1:knon)= run_off_lic(1:knon) - call gath2cpl(bidule, tmp_rlic, klon, knon, iim, jjm, knindex) - - ! calcul albedo - - CALL albsno(klon, knon, dtime, agesno, alb_neig, precip_snow) - WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0. - zfra(1:knon) = MAX(0.0, MIN(1.0, snow(1:knon)/(snow(1:knon)+10.0))) - alb_new(1 : knon) = alb_neig(1 : knon)*zfra(1:knon) + & - 0.6 * (1.0-zfra(1:knon)) - - !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux" - ! alb_new(1 : knon) = 0.6 !IM cf FH/GK - ! alb_new(1 : knon) = 0.82 - ! alb_new(1 : knon) = 0.77 !211003 Ksta0.77 - ! alb_new(1 : knon) = 0.8 !KstaTER0.8 & LMD_ARMIP5 - !IM: KstaTER0.77 & LMD_ARMIP6 - alb_new(1 : knon) = 0.77 - - - ! Rugosite - - z0_new = rugoro - - ! Remplissage des pourcentages de surface - - pctsrf_new(:, nisurf) = pctsrf(:, nisurf) - - alblw(1:knon) = alb_new(1:knon) - else - write(*, *)'Index surface = ', nisurf - abort_message = 'Index surface non valable' - call abort_gcm(modname, abort_message, 1) - endif - - END SUBROUTINE interfsurf_hq - - !************************ - - SUBROUTINE interfoce_slab(klon, debut, itap, dtime, ijour, & - radsol, fluxo, fluxg, pctsrf, & - tslab, seaice, pctsrf_slab) - - ! Cette routine calcule la temperature d'un slab ocean, la glace de mer - ! et les pourcentages de la maille couverte par l'ocean libre et/ou - ! la glace de mer pour un "slab" ocean de 50m - - ! I. Musat 04.02.2005 - - ! input: - ! klon nombre total de points de grille - ! debut logical: 1er appel a la physique - ! itap numero du pas de temps - ! dtime pas de temps de la physique (en s) - ! ijour jour dans l'annee en cours - ! radsol rayonnement net au sol (LW + SW) - ! fluxo flux turbulent (sensible + latent) sur les mailles oceaniques - ! fluxg flux de conduction entre la surface de la glace de mer et l'ocean - ! pctsrf tableau des pourcentages de surface de chaque maille - ! output: - ! tslab temperature de l'ocean libre - ! seaice glace de mer (kg/m2) - ! pctsrf_slab "pourcentages" (valeurs entre 0. et 1.) surfaces issus du slab - - use indicesol - use clesphys - use abort_gcm_m, only: abort_gcm - use SUPHEC_M - - ! Parametres d'entree - integer, intent(IN) :: klon - logical, intent(IN) :: debut - INTEGER, intent(IN) :: itap - REAL, intent(IN) :: dtime - INTEGER, intent(IN) :: ijour - REAL, dimension(klon), intent(IN) :: radsol - REAL, dimension(klon), intent(IN) :: fluxo - REAL, dimension(klon), intent(IN) :: fluxg - real, dimension(klon, nbsrf), intent(IN) :: pctsrf - ! Parametres de sortie - real, dimension(klon), intent(INOUT) :: tslab - real, dimension(klon), intent(INOUT) :: seaice ! glace de mer (kg/m2) - real, dimension(klon, nbsrf), intent(OUT) :: pctsrf_slab - - ! Variables locales : - INTEGER, save :: lmt_pas, julien, idayvrai - REAL, parameter :: unjour=86400. - real, allocatable, dimension(:), save :: tmp_tslab, tmp_seaice - REAL, allocatable, dimension(:), save :: slab_bils - REAL, allocatable, dimension(:), save :: lmt_bils - logical, save :: check = .false. - - REAL, parameter :: cyang=50.0 * 4.228e+06 ! capacite calorifique volumetrique de l'eau J/(m2 K) - REAL, parameter :: cbing=0.334e+05 ! J/kg - real, dimension(klon) :: siceh !hauteur de la glace de mer (m) - INTEGER :: i - integer :: sum_error, error - REAL :: zz, za, zb - - character (len = 80) :: abort_message - character (len = 20) :: modname = 'interfoce_slab' - - julien = MOD(ijour, 360) - sum_error = 0 - IF (debut) THEN - allocate(slab_bils(klon), stat = error); sum_error = sum_error + error - allocate(lmt_bils(klon), stat = error); sum_error = sum_error + error - allocate(tmp_tslab(klon), stat = error); sum_error = sum_error + error - allocate(tmp_seaice(klon), stat = error); sum_error = sum_error + error - if (sum_error /= 0) then - abort_message='Pb allocation var. slab_bils, lmt_bils, tmp_tslab, tmp_seaice' - call abort_gcm(modname, abort_message, 1) - endif - tmp_tslab=tslab - tmp_seaice=seaice - lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour - - IF (check) THEN - PRINT*, 'interfoce_slab klon, debut, itap, dtime, ijour, & - & lmt_pas ', klon, debut, itap, dtime, ijour, & - lmt_pas - ENDIF !check - - PRINT*, '************************' - PRINT*, 'SLAB OCEAN est actif, prenez precautions !' - PRINT*, '************************' - - ! a mettre un slab_bils aussi en force !!! - - DO i = 1, klon - slab_bils(i) = 0.0 - ENDDO - - ENDIF !debut - pctsrf_slab(1:klon, 1:nbsrf) = pctsrf(1:klon, 1:nbsrf) - - ! lecture du bilan au sol lmt_bils issu d'une simulation forcee en debut de journee - - IF (MOD(itap, lmt_pas) .EQ. 1) THEN !1er pas de temps de la journee - idayvrai = ijour - CALL condsurf(julien, idayvrai, lmt_bils) - ENDIF !(MOD(itap-1, lmt_pas) .EQ. 0) THEN - - DO i = 1, klon - IF((pctsrf_slab(i, is_oce).GT.epsfra).OR. & - (pctsrf_slab(i, is_sic).GT.epsfra)) THEN - - ! fabriquer de la glace si congelation atteinte: - - IF (tmp_tslab(i).LT.(RTT-1.8)) THEN - zz = (RTT-1.8)-tmp_tslab(i) - tmp_seaice(i) = tmp_seaice(i) + cyang/cbing * zz - seaice(i) = tmp_seaice(i) - tmp_tslab(i) = RTT-1.8 - ENDIF - - ! faire fondre de la glace si temperature est superieure a 0: - - IF ((tmp_tslab(i).GT.RTT) .AND. (tmp_seaice(i).GT.0.0)) THEN - zz = cyang/cbing * (tmp_tslab(i)-RTT) - zz = MIN(zz, tmp_seaice(i)) - tmp_seaice(i) = tmp_seaice(i) - zz - seaice(i) = tmp_seaice(i) - tmp_tslab(i) = tmp_tslab(i) - zz*cbing/cyang - ENDIF - - ! limiter la glace de mer a 10 metres (10000 kg/m2) - - IF(tmp_seaice(i).GT.45.) THEN - tmp_seaice(i) = MIN(tmp_seaice(i), 10000.0) - ELSE - tmp_seaice(i) = 0. - ENDIF - seaice(i) = tmp_seaice(i) - siceh(i)=tmp_seaice(i)/1000. !en metres - - ! determiner la nature du sol (glace de mer ou ocean libre): - - ! on fait dependre la fraction de seaice "pctsrf(i, is_sic)" - ! de l'epaisseur de seaice : - ! pctsrf(i, is_sic)=1. si l'epaisseur de la glace de mer est >= a 20cm - ! et pctsrf(i, is_sic) croit lineairement avec seaice de 0. a 20cm d'epaisseur - - pctsrf_slab(i, is_sic)=MIN(siceh(i)/0.20, & - 1.-(pctsrf_slab(i, is_ter)+pctsrf_slab(i, is_lic))) - pctsrf_slab(i, is_oce)=1.0 - & - (pctsrf_slab(i, is_ter)+pctsrf_slab(i, is_lic)+pctsrf_slab(i, is_sic)) - ENDIF !pctsrf - ENDDO - - ! Calculer le bilan du flux de chaleur au sol : - - DO i = 1, klon - za = radsol(i) + fluxo(i) - zb = fluxg(i) - IF((pctsrf_slab(i, is_oce).GT.epsfra).OR. & - (pctsrf_slab(i, is_sic).GT.epsfra)) THEN - slab_bils(i)=slab_bils(i)+(za*pctsrf_slab(i, is_oce) & - +zb*pctsrf_slab(i, is_sic))/ FLOAT(lmt_pas) - ENDIF - ENDDO !klon - - ! calcul tslab - - IF (MOD(itap, lmt_pas).EQ.0) THEN !fin de journee - DO i = 1, klon - IF ((pctsrf_slab(i, is_oce).GT.epsfra).OR. & - (pctsrf_slab(i, is_sic).GT.epsfra)) THEN - tmp_tslab(i) = tmp_tslab(i) + & - (slab_bils(i)-lmt_bils(i)) & - /cyang*unjour - ! on remet l'accumulation a 0 - slab_bils(i) = 0. - ENDIF !pctsrf - ENDDO !klon - ENDIF !(MOD(itap, lmt_pas).EQ.0) THEN - - tslab = tmp_tslab - seaice = tmp_seaice - END SUBROUTINE interfoce_slab - - !************************ - - SUBROUTINE interfoce_lim(itime, dtime, jour, & - klon, nisurf, knon, knindex, & - debut, & - lmt_sst, pctsrf_new) - - ! Cette routine sert d'interface entre le modele atmospherique et - ! un fichier de conditions aux limites - - ! L. Fairhead 02/2000 - - use abort_gcm_m, only: abort_gcm - use indicesol - - integer, intent(IN) :: itime ! numero du pas de temps courant - real , intent(IN) :: dtime ! pas de temps de la physique (en s) - integer, intent(IN) :: jour ! jour a lire dans l'annee - integer, intent(IN) :: nisurf ! index de la surface a traiter (1 = sol continental) - integer, intent(IN) :: knon ! nombre de points dans le domaine a traiter - integer, intent(IN) :: klon ! taille de la grille - integer, dimension(klon), intent(in) :: knindex ! index des points de la surface a traiter - logical, intent(IN) :: debut ! logical: 1er appel a la physique (initialisation) - - ! Parametres de sortie - ! output: - ! lmt_sst SST lues dans le fichier de CL - ! pctsrf_new sous-maille fractionnelle - real, intent(out), dimension(klon) :: lmt_sst - real, intent(out), dimension(klon, nbsrf) :: pctsrf_new - - ! Variables locales - integer :: ii - INTEGER, save :: lmt_pas ! frequence de lecture des conditions limites - ! (en pas de physique) - logical, save :: deja_lu ! pour indiquer que le jour a lire a deja - ! lu pour une surface precedente - integer, save :: jour_lu - integer :: ierr - character (len = 20) :: modname = 'interfoce_lim' - character (len = 80) :: abort_message - logical, save :: newlmt = .TRUE. - logical, save :: check = .FALSE. - ! Champs lus dans le fichier de CL - real, allocatable , save, dimension(:) :: sst_lu, rug_lu, nat_lu - real, allocatable , save, dimension(:, :) :: pct_tmp - - ! quelques variables pour netcdf - - include "netcdf.inc" - integer :: nid, nvarid - integer, dimension(2) :: start, epais - - ! -------------------------------------------------- - - if (debut .and. .not. allocated(sst_lu)) then - lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour - jour_lu = jour - 1 - allocate(sst_lu(klon)) - allocate(nat_lu(klon)) - allocate(pct_tmp(klon, nbsrf)) - endif - - if ((jour - jour_lu) /= 0) deja_lu = .false. - - if (check) write(*, *)modname, ' :: jour, jour_lu, deja_lu', jour, jour_lu, & - deja_lu - if (check) write(*, *)modname, ' :: itime, lmt_pas ', itime, lmt_pas, dtime - - ! Tester d'abord si c'est le moment de lire le fichier - if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu) then - - ! Ouverture du fichier - - ierr = NF_OPEN ('limit.nc', NF_NOWRITE, nid) - if (ierr.NE.NF_NOERR) then - abort_message & - = 'Pb d''ouverture du fichier de conditions aux limites' - call abort_gcm(modname, abort_message, 1) - endif - - ! La tranche de donnees a lire: - - start(1) = 1 - start(2) = jour - epais(1) = klon - epais(2) = 1 - - if (newlmt) then - - ! Fraction "ocean" - - ierr = NF_INQ_VARID(nid, 'FOCE', nvarid) - if (ierr /= NF_NOERR) then - abort_message = 'Le champ est absent' - call abort_gcm(modname, abort_message, 1) - endif - ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, pct_tmp(1, is_oce)) - if (ierr /= NF_NOERR) then - abort_message = 'Lecture echouee pour ' - call abort_gcm(modname, abort_message, 1) - endif - - ! Fraction "glace de mer" - - ierr = NF_INQ_VARID(nid, 'FSIC', nvarid) - if (ierr /= NF_NOERR) then - abort_message = 'Le champ est absent' - call abort_gcm(modname, abort_message, 1) - endif - ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, pct_tmp(1, is_sic)) - if (ierr /= NF_NOERR) then - abort_message = 'Lecture echouee pour ' - call abort_gcm(modname, abort_message, 1) - endif - - ! Fraction "terre" - - ierr = NF_INQ_VARID(nid, 'FTER', nvarid) - if (ierr /= NF_NOERR) then - abort_message = 'Le champ est absent' - call abort_gcm(modname, abort_message, 1) - endif - ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, pct_tmp(1, is_ter)) - if (ierr /= NF_NOERR) then - abort_message = 'Lecture echouee pour ' - call abort_gcm(modname, abort_message, 1) - endif - - ! Fraction "glacier terre" - - ierr = NF_INQ_VARID(nid, 'FLIC', nvarid) - if (ierr /= NF_NOERR) then - abort_message = 'Le champ est absent' - call abort_gcm(modname, abort_message, 1) - endif - ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, pct_tmp(1, is_lic)) - if (ierr /= NF_NOERR) then - abort_message = 'Lecture echouee pour ' - call abort_gcm(modname, abort_message, 1) - endif - - else ! on en est toujours a rnatur - - ierr = NF_INQ_VARID(nid, 'NAT', nvarid) - if (ierr /= NF_NOERR) then - abort_message = 'Le champ est absent' - call abort_gcm(modname, abort_message, 1) - endif - ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, nat_lu) - if (ierr /= NF_NOERR) then - abort_message = 'Lecture echouee pour ' - call abort_gcm(modname, abort_message, 1) - endif - - ! Remplissage des fractions de surface - ! nat = 0, 1, 2, 3 pour ocean, terre, glacier, seaice - - pct_tmp = 0.0 - do ii = 1, klon - pct_tmp(ii, nint(nat_lu(ii)) + 1) = 1. - enddo - - - ! On se retrouve avec ocean en 1 et terre en 2 alors qu'on veut le contraire - - pctsrf_new = pct_tmp - pctsrf_new (:, 2)= pct_tmp (:, 1) - pctsrf_new (:, 1)= pct_tmp (:, 2) - pct_tmp = pctsrf_new - endif ! fin test sur newlmt - - ! Lecture SST - - ierr = NF_INQ_VARID(nid, 'SST', nvarid) - if (ierr /= NF_NOERR) then - abort_message = 'Le champ est absent' - call abort_gcm(modname, abort_message, 1) - endif - ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, sst_lu) - if (ierr /= NF_NOERR) then - abort_message = 'Lecture echouee pour ' - call abort_gcm(modname, abort_message, 1) - endif - - - ! Fin de lecture - - ierr = NF_CLOSE(nid) - deja_lu = .true. - jour_lu = jour - endif - - ! Recopie des variables dans les champs de sortie - - lmt_sst = 999999999. - do ii = 1, knon - lmt_sst(ii) = sst_lu(knindex(ii)) - enddo - - pctsrf_new(:, is_oce) = pct_tmp(:, is_oce) - pctsrf_new(:, is_sic) = pct_tmp(:, is_sic) - - END SUBROUTINE interfoce_lim - - !************************ - - SUBROUTINE interfsur_lim(itime, dtime, jour, & - klon, nisurf, knon, knindex, & - debut, & - lmt_alb, lmt_rug) - - ! Cette routine sert d'interface entre le modèle atmosphérique et - ! un fichier de conditions aux limites. - - ! L. Fairhead 02/2000 - - use abort_gcm_m, only: abort_gcm - - ! Parametres d'entree - ! input: - ! itime numero du pas de temps courant - ! dtime pas de temps de la physique (en s) - ! jour jour a lire dans l'annee - ! nisurf index de la surface a traiter (1 = sol continental) - ! knon nombre de points dans le domaine a traiter - ! knindex index des points de la surface a traiter - ! klon taille de la grille - ! debut logical: 1er appel a la physique (initialisation) - integer, intent(IN) :: itime - real , intent(IN) :: dtime - integer, intent(IN) :: jour - integer, intent(IN) :: nisurf - integer, intent(IN) :: knon - integer, intent(IN) :: klon - integer, dimension(klon), intent(in) :: knindex - logical, intent(IN) :: debut - - ! Parametres de sortie - ! output: - ! lmt_sst SST lues dans le fichier de CL - ! lmt_alb Albedo lu - ! lmt_rug longueur de rugosité lue - ! pctsrf_new sous-maille fractionnelle - real, intent(out), dimension(klon) :: lmt_alb - real, intent(out), dimension(klon) :: lmt_rug - - ! Variables locales - integer :: ii - integer, save :: lmt_pas ! frequence de lecture des conditions limites - ! (en pas de physique) - logical, save :: deja_lu_sur! pour indiquer que le jour a lire a deja - ! lu pour une surface precedente - integer, save :: jour_lu_sur - integer :: ierr - character (len = 20) :: modname = 'interfsur_lim' - character (len = 80) :: abort_message - logical, save :: newlmt = .false. - logical, save :: check = .false. - ! Champs lus dans le fichier de CL - real, allocatable , save, dimension(:) :: alb_lu, rug_lu - - ! quelques variables pour netcdf - - include "netcdf.inc" - integer , save :: nid, nvarid - integer, dimension(2), save :: start, epais - - !------------------------------------------------------------ - - if (debut) then - lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour - jour_lu_sur = jour - 1 - allocate(alb_lu(klon)) - allocate(rug_lu(klon)) - endif - - if ((jour - jour_lu_sur) /= 0) deja_lu_sur = .false. - - if (check) write(*, *)modname, ':: jour_lu_sur, deja_lu_sur', jour_lu_sur, & - deja_lu_sur - if (check) write(*, *)modname, ':: itime, lmt_pas', itime, lmt_pas - - ! Tester d'abord si c'est le moment de lire le fichier - if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu_sur) then - - ! Ouverture du fichier - - ierr = NF_OPEN ('limit.nc', NF_NOWRITE, nid) - if (ierr.NE.NF_NOERR) then - abort_message & - = 'Pb d''ouverture du fichier de conditions aux limites' - call abort_gcm(modname, abort_message, 1) - endif - - ! La tranche de donnees a lire: - - start(1) = 1 - start(2) = jour - epais(1) = klon - epais(2) = 1 - - ! Lecture Albedo - - ierr = NF_INQ_VARID(nid, 'ALB', nvarid) - if (ierr /= NF_NOERR) then - abort_message = 'Le champ est absent' - call abort_gcm(modname, abort_message, 1) - endif - ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, alb_lu) - if (ierr /= NF_NOERR) then - abort_message = 'Lecture echouee pour ' - call abort_gcm(modname, abort_message, 1) - endif - - ! Lecture rugosité - - ierr = NF_INQ_VARID(nid, 'RUG', nvarid) - if (ierr /= NF_NOERR) then - abort_message = 'Le champ est absent' - call abort_gcm(modname, abort_message, 1) - endif - ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, rug_lu) - if (ierr /= NF_NOERR) then - abort_message = 'Lecture echouee pour ' - call abort_gcm(modname, abort_message, 1) - endif - - - ! Fin de lecture - - ierr = NF_CLOSE(nid) - deja_lu_sur = .true. - jour_lu_sur = jour - endif - - ! Recopie des variables dans les champs de sortie - -!!$ lmt_alb = 0.0 -!!$ lmt_rug = 0.0 - lmt_alb = 999999. - lmt_rug = 999999. - DO ii = 1, knon - lmt_alb(ii) = alb_lu(knindex(ii)) - lmt_rug(ii) = rug_lu(knindex(ii)) - enddo - - END SUBROUTINE interfsur_lim - - !************************ - - SUBROUTINE calcul_fluxs( klon, knon, nisurf, dtime, & - tsurf, p1lay, cal, beta, coef1lay, ps, & - precip_rain, precip_snow, snow, qsurf, & - radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, & - petAcoef, peqAcoef, petBcoef, peqBcoef, & - tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) - - ! Cette routine calcule les fluxs en h et q a l'interface et eventuellement - ! une temperature de surface (au cas ou ok_veget = false) - - ! L. Fairhead 4/2000 - - ! input: - ! knon nombre de points a traiter - ! nisurf surface a traiter - ! tsurf temperature de surface - ! p1lay pression 1er niveau (milieu de couche) - ! cal capacite calorifique du sol - ! beta evap reelle - ! coef1lay coefficient d'echange - ! ps pression au sol - ! precip_rain precipitations liquides - ! precip_snow precipitations solides - ! snow champs hauteur de neige - ! runoff runoff en cas de trop plein - ! petAcoef coeff. A de la resolution de la CL pour t - ! peqAcoef coeff. A de la resolution de la CL pour q - ! petBcoef coeff. B de la resolution de la CL pour t - ! peqBcoef coeff. B de la resolution de la CL pour q - ! radsol rayonnement net aus sol (LW + SW) - ! dif_grnd coeff. diffusion vers le sol profond - - ! output: - ! tsurf_new temperature au sol - ! qsurf humidite de l'air au dessus du sol - ! fluxsens flux de chaleur sensible - ! fluxlat flux de chaleur latente - ! dflux_s derivee du flux de chaleur sensible / Ts - ! dflux_l derivee du flux de chaleur latente / Ts - - - use indicesol - use abort_gcm_m, only: abort_gcm - use yoethf_m - use fcttre, only: thermcep, foeew, qsats, qsatl, foede, dqsats, dqsatl - use SUPHEC_M - - ! Parametres d'entree - integer, intent(IN) :: knon, nisurf, klon - real , intent(IN) :: dtime - real, dimension(klon), intent(IN) :: petAcoef, peqAcoef - real, dimension(klon), intent(IN) :: petBcoef, peqBcoef - real, dimension(klon), intent(IN) :: ps, q1lay - real, dimension(klon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay - real, dimension(klon), intent(IN) :: precip_rain, precip_snow - real, dimension(klon), intent(IN) :: radsol, dif_grnd - real, dimension(klon), intent(IN) :: t1lay, u1lay, v1lay - real, dimension(klon), intent(INOUT) :: snow, qsurf - - ! Parametres sorties - real, dimension(klon), intent(OUT):: tsurf_new, evap, fluxsens, fluxlat - real, dimension(klon), intent(OUT):: dflux_s, dflux_l - - ! Variables locales - integer :: i - real, dimension(klon) :: zx_mh, zx_nh, zx_oh - real, dimension(klon) :: zx_mq, zx_nq, zx_oq - real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef - real, dimension(klon) :: zx_sl, zx_k1 - real, dimension(klon) :: zx_q_0 , d_ts - real :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh - real :: bilan_f, fq_fonte - REAL :: subli, fsno - REAL :: qsat_new, q1_new - real, parameter :: t_grnd = 271.35, t_coup = 273.15 - !! PB temporaire en attendant mieux pour le modele de neige - REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15) - - logical, save :: check = .false. - character (len = 20) :: modname = 'calcul_fluxs' - logical, save :: fonte_neige = .false. - real, save :: max_eau_sol = 150.0 - character (len = 80) :: abort_message - logical, save :: first = .true., second=.false. - - if (check) write(*, *)'Entree ', modname, ' surface = ', nisurf - - IF (check) THEN - WRITE(*, *)' radsol (min, max)' & - , MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon)) - !!CALL flush(6) - ENDIF - - if (size(coastalflow) /= knon .AND. nisurf == is_ter) then - write(*, *)'Bizarre, le nombre de points continentaux' - write(*, *)'a change entre deux appels. J''arrete ...' - abort_message='Pb run_off' - call abort_gcm(modname, abort_message, 1) - endif - - ! Traitement neige et humidite du sol - - ! Initialisation - - evap = 0. - fluxsens=0. - fluxlat=0. - dflux_s = 0. - dflux_l = 0. - - ! zx_qs = qsat en kg/kg - - DO i = 1, knon - zx_pkh(i) = (ps(i)/ps(i))**RKAPPA - IF (thermcep) THEN - zdelta=MAX(0., SIGN(1., rtt-tsurf(i))) - zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta - zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i)) - zx_qs= r2es * FOEEW(tsurf(i), zdelta)/ps(i) - zx_qs=MIN(0.5, zx_qs) - zcor=1./(1.-retv*zx_qs) - zx_qs=zx_qs*zcor - zx_dq_s_dh = FOEDE(tsurf(i), zdelta, zcvm5, zx_qs, zcor) & - /RLVTT / zx_pkh(i) - ELSE - IF (tsurf(i).LT.t_coup) THEN - zx_qs = qsats(tsurf(i)) / ps(i) - zx_dq_s_dh = dqsats(tsurf(i), zx_qs)/RLVTT & - / zx_pkh(i) - ELSE - zx_qs = qsatl(tsurf(i)) / ps(i) - zx_dq_s_dh = dqsatl(tsurf(i), zx_qs)/RLVTT & - / zx_pkh(i) - ENDIF - ENDIF - zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh - zx_qsat(i) = zx_qs - zx_coef(i) = coef1lay(i) & - * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) & - * p1lay(i)/(RD*t1lay(i)) - - ENDDO - - ! === Calcul de la temperature de surface === - - ! zx_sl = chaleur latente d'evaporation ou de sublimation - - do i = 1, knon - zx_sl(i) = RLVTT - if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT - zx_k1(i) = zx_coef(i) - enddo - - do i = 1, knon - ! Q - zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime) - zx_mq(i) = beta(i) * zx_k1(i) * & - (peqAcoef(i) - zx_qsat(i) & - + zx_dq_s_dt(i) * tsurf(i)) & - / zx_oq(i) - zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) & - / zx_oq(i) - - ! H - zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime) - zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i) - zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i) - - ! Tsurface - tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * & - (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) & - + dif_grnd(i) * t_grnd * dtime)/ & - ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( & - zx_nh(i) + zx_sl(i) * zx_nq(i)) & - + dtime * dif_grnd(i)) - - - ! Y'a-t-il fonte de neige? - - ! fonte_neige = (nisurf /= is_oce) .AND. & - ! & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) & - ! & .AND. (tsurf_new(i) >= RTT) - ! if (fonte_neige) tsurf_new(i) = RTT - d_ts(i) = tsurf_new(i) - tsurf(i) - ! zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i) - ! zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i) - !== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas - !== flux_t est le flux de cpt (energie sensible): j/(m**2 s) - evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i) - fluxlat(i) = - evap(i) * zx_sl(i) - fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i) - ! Derives des flux dF/dTs (W m-2 K-1): - dflux_s(i) = zx_nh(i) - dflux_l(i) = (zx_sl(i) * zx_nq(i)) - ! Nouvelle valeure de l'humidite au dessus du sol - qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i) - q1_new = peqAcoef(i) - peqBcoef(i)*evap(i)*dtime - qsurf(i)=q1_new*(1.-beta(i)) + beta(i)*qsat_new - ENDDO - - END SUBROUTINE calcul_fluxs - - !************************ - - SUBROUTINE fonte_neige( klon, knon, nisurf, dtime, & - tsurf, p1lay, cal, beta, coef1lay, ps, & - precip_rain, precip_snow, snow, qsol, & - radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, & - petAcoef, peqAcoef, petBcoef, peqBcoef, & - tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & - fqcalving, ffonte, run_off_lic_0) - - ! Routine de traitement de la fonte de la neige dans le cas du traitement - ! de sol simplifié - - ! LF 03/2001 - ! input: - ! knon nombre de points a traiter - ! nisurf surface a traiter - ! tsurf temperature de surface - ! p1lay pression 1er niveau (milieu de couche) - ! cal capacite calorifique du sol - ! beta evap reelle - ! coef1lay coefficient d'echange - ! ps pression au sol - ! precip_rain precipitations liquides - ! precip_snow precipitations solides - ! snow champs hauteur de neige - ! qsol hauteur d'eau contenu dans le sol - ! runoff runoff en cas de trop plein - ! petAcoef coeff. A de la resolution de la CL pour t - ! peqAcoef coeff. A de la resolution de la CL pour q - ! petBcoef coeff. B de la resolution de la CL pour t - ! peqBcoef coeff. B de la resolution de la CL pour q - ! radsol rayonnement net aus sol (LW + SW) - ! dif_grnd coeff. diffusion vers le sol profond - - ! output: - ! tsurf_new temperature au sol - ! fluxsens flux de chaleur sensible - ! fluxlat flux de chaleur latente - ! dflux_s derivee du flux de chaleur sensible / Ts - ! dflux_l derivee du flux de chaleur latente / Ts - ! in/out: - ! run_off_lic_0 run off glacier du pas de temps précedent - - - use indicesol - use SUPHEC_M - use yoethf_m - use fcttre - !IM cf JLD - - ! Parametres d'entree - integer, intent(IN) :: knon, nisurf, klon - real , intent(IN) :: dtime - real, dimension(klon), intent(IN) :: petAcoef, peqAcoef - real, dimension(klon), intent(IN) :: petBcoef, peqBcoef - real, dimension(klon), intent(IN) :: ps, q1lay - real, dimension(klon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay - real, dimension(klon), intent(IN) :: precip_rain, precip_snow - real, dimension(klon), intent(IN) :: radsol, dif_grnd - real, dimension(klon), intent(IN) :: t1lay, u1lay, v1lay - real, dimension(klon), intent(INOUT) :: snow, qsol - - ! Parametres sorties - real, dimension(klon), intent(INOUT):: tsurf_new, evap, fluxsens, fluxlat - real, dimension(klon), intent(INOUT):: dflux_s, dflux_l - ! Flux thermique utiliser pour fondre la neige - real, dimension(klon), intent(INOUT):: ffonte - ! Flux d'eau "perdue" par la surface et necessaire pour que limiter la - ! hauteur de neige, en kg/m2/s - real, dimension(klon), intent(INOUT):: fqcalving - real, dimension(klon), intent(INOUT):: run_off_lic_0 - ! Variables locales - ! Masse maximum de neige (kg/m2). Au dessus de ce seuil, la neige - ! en exces "s'ecoule" (calving) - ! real, parameter :: snow_max=1. - !IM cf JLD/GK - real, parameter :: snow_max=3000. - integer :: i - real, dimension(klon) :: zx_mh, zx_nh, zx_oh - real, dimension(klon) :: zx_mq, zx_nq, zx_oq - real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef - real, dimension(klon) :: zx_sl, zx_k1 - real, dimension(klon) :: zx_q_0 , d_ts - real :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh - real :: bilan_f, fq_fonte - REAL :: subli, fsno - REAL, DIMENSION(klon) :: bil_eau_s, snow_evap - real, parameter :: t_grnd = 271.35, t_coup = 273.15 - !! PB temporaire en attendant mieux pour le modele de neige - ! REAL, parameter :: chasno = RLMLT/(2.3867E+06*0.15) - REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15) - !IM cf JLD/ GKtest - REAL, parameter :: chaice = 3.334E+05/(2.3867E+06*0.15) - ! fin GKtest - - logical, save :: check = .FALSE. - character (len = 20) :: modname = 'fonte_neige' - logical, save :: neige_fond = .false. - real, save :: max_eau_sol = 150.0 - character (len = 80) :: abort_message - logical, save :: first = .true., second=.false. - real :: coeff_rel - - if (check) write(*, *)'Entree ', modname, ' surface = ', nisurf - - ! Initialisations - coeff_rel = dtime/(tau_calv * rday) - bil_eau_s = 0. - DO i = 1, knon - zx_pkh(i) = (ps(i)/ps(i))**RKAPPA - IF (thermcep) THEN - zdelta=MAX(0., SIGN(1., rtt-tsurf(i))) - zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta - zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i)) - zx_qs= r2es * FOEEW(tsurf(i), zdelta)/ps(i) - zx_qs=MIN(0.5, zx_qs) - zcor=1./(1.-retv*zx_qs) - zx_qs=zx_qs*zcor - zx_dq_s_dh = FOEDE(tsurf(i), zdelta, zcvm5, zx_qs, zcor) & - /RLVTT / zx_pkh(i) - ELSE - IF (tsurf(i).LT.t_coup) THEN - zx_qs = qsats(tsurf(i)) / ps(i) - zx_dq_s_dh = dqsats(tsurf(i), zx_qs)/RLVTT & - / zx_pkh(i) - ELSE - zx_qs = qsatl(tsurf(i)) / ps(i) - zx_dq_s_dh = dqsatl(tsurf(i), zx_qs)/RLVTT & - / zx_pkh(i) - ENDIF - ENDIF - zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh - zx_qsat(i) = zx_qs - zx_coef(i) = coef1lay(i) & - * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) & - * p1lay(i)/(RD*t1lay(i)) - ENDDO - - ! === Calcul de la temperature de surface === - - ! zx_sl = chaleur latente d'evaporation ou de sublimation - - do i = 1, knon - zx_sl(i) = RLVTT - if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT - zx_k1(i) = zx_coef(i) - enddo - - do i = 1, knon - ! Q - zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime) - zx_mq(i) = beta(i) * zx_k1(i) * & - (peqAcoef(i) - zx_qsat(i) & - + zx_dq_s_dt(i) * tsurf(i)) & - / zx_oq(i) - zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) & - / zx_oq(i) - - ! H - zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime) - zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i) - zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i) - enddo - - WHERE (precip_snow > 0.) snow = snow + (precip_snow * dtime) - snow_evap = 0. - WHERE (evap > 0. ) - snow_evap = MIN (snow / dtime, evap) - snow = snow - snow_evap * dtime - snow = MAX(0.0, snow) - end where - - ! bil_eau_s = bil_eau_s + (precip_rain * dtime) - (evap - snow_evap) * dtime - bil_eau_s = (precip_rain * dtime) - (evap - snow_evap) * dtime - - - ! Y'a-t-il fonte de neige? - - ffonte=0. - do i = 1, knon - neige_fond = ((snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) & - .AND. tsurf_new(i) >= RTT) - if (neige_fond) then - fq_fonte = MIN( MAX((tsurf_new(i)-RTT )/chasno, 0.0), snow(i)) - ffonte(i) = fq_fonte * RLMLT/dtime - snow(i) = max(0., snow(i) - fq_fonte) - bil_eau_s(i) = bil_eau_s(i) + fq_fonte - tsurf_new(i) = tsurf_new(i) - fq_fonte * chasno - !IM cf JLD OK - !IM cf JLD/ GKtest fonte aussi pour la glace - IF (nisurf == is_sic .OR. nisurf == is_lic ) THEN - fq_fonte = MAX((tsurf_new(i)-RTT )/chaice, 0.0) - ffonte(i) = ffonte(i) + fq_fonte * RLMLT/dtime - bil_eau_s(i) = bil_eau_s(i) + fq_fonte - tsurf_new(i) = RTT - ENDIF - d_ts(i) = tsurf_new(i) - tsurf(i) - endif - - ! s'il y a une hauteur trop importante de neige, elle s'coule - fqcalving(i) = max(0., snow(i) - snow_max)/dtime - snow(i)=min(snow(i), snow_max) - - IF (nisurf == is_ter) then - qsol(i) = qsol(i) + bil_eau_s(i) - run_off(i) = run_off(i) + MAX(qsol(i) - max_eau_sol, 0.0) - qsol(i) = MIN(qsol(i), max_eau_sol) - else if (nisurf == is_lic) then - run_off_lic(i) = (coeff_rel * fqcalving(i)) + & - (1. - coeff_rel) * run_off_lic_0(i) - run_off_lic_0(i) = run_off_lic(i) - run_off_lic(i) = run_off_lic(i) + bil_eau_s(i)/dtime - endif - enddo - - END SUBROUTINE fonte_neige - END MODULE interface_surf