subroutine hydrol(ngrid,nq,ptimestep,rnat,tsurf, & qsurf,dqsurf,dqs_hyd,pcapcal, & albedo,albedo_bareground, & albedo_snow_SPECTV,albedo_co2_ice_SPECTV, & mu0,pdtsurf,pdtsurf_hyd,hice, & pctsrf_sic,sea_ice) use ioipsl_getin_p_mod, only: getin_p use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid use mod_phys_lmdz_para, only : is_master, gather, scatter use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol USE surfdat_h use comdiurn_h USE geometry_mod, only: cell_area USE tracer_h ! use slab_ice_h USE ocean_slab_mod, ONLY: alb_ice_min,h_alb_ice,snow_min use callkeys_mod, only: albedosnow,alb_ocean,albedoco2ice,ok_slab_ocean,Tsaldiff,maxicethick,co2cond use radinc_h, only : L_NSPECTV implicit none !================================================================== ! ! Purpose ! ------- ! Calculate the surface hydrology and albedo changes. ! Both for oceanic and continental regions ! ! Authors ! ------- ! Adapted from LMDTERRE by B. Charnay (2010). Further ! Modifications by R. Wordsworth (2010). ! Spectral albedo by M. Turbet (2015). ! ! Called by ! --------- ! physiq.F ! ! Calls ! ----- ! none ! ! Notes ! ----- ! rnat is terrain type: 0-ocean; 1-continent ! !================================================================== integer ngrid,nq ! Inputs ! ------ real snowlayer parameter (snowlayer=33.0) ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm real oceantime ! this is a relaxation timescale for the oceanbulkavg parameterization parameter (oceantime=10*24*3600) logical,save :: oceanbulkavg ! simple parameterization to relax ocean temperatures to the global mean value (crude, 0th order parameterization to mimick ocean heat transport) logical,save :: activerunoff ! enable simple runoff scheme? logical,save :: oceanalbvary ! simple parameterization to account for the effect of solar zenith angle on the ocean albedo (for the moment it is not used, but to be included in the future) !$OMP THREADPRIVATE(oceanbulkavg,activerunoff,oceanalbvary) ! Arguments ! --------- real rnat(ngrid) ! rnat is terrain type: 0-ocean; 1-continent real,dimension(:),allocatable,save :: runoff real totalrunoff, tsea, oceanarea save oceanarea !$OMP THREADPRIVATE(runoff,oceanarea) real ptimestep real mu0(ngrid) real qsurf(ngrid,nq), tsurf(ngrid) real dqsurf(ngrid,nq), pdtsurf(ngrid) real hice(ngrid) real albedo(ngrid,L_NSPECTV) real albedo_bareground(ngrid) real albedo_snow_SPECTV(L_NSPECTV) real albedo_co2_ice_SPECTV(L_NSPECTV) real pctsrf_sic(ngrid), sea_ice(ngrid) real oceanarea2 ! Output ! ------ real dqs_hyd(ngrid,nq) real pdtsurf_hyd(ngrid) ! Local ! ----- real a,b,E integer ig,iq, nw real fsnoi, subli, fauxo real twater(ngrid) real pcapcal(ngrid) real hicebis(ngrid) real zqsurf(ngrid,nq) real ztsurf(ngrid) real albedo_sic, alb_ice real frac_snow integer, save :: ivap, iliq, iice !$OMP THREADPRIVATE(ivap,iliq,iice) logical, save :: firstcall=.true. !$OMP THREADPRIVATE(firstcall) real :: runoffamount(ngrid) !#ifdef CPP_PARA real :: runoffamount_glo(klon_glo) real :: zqsurf_iliq_glo(klon_glo) real :: rnat_glo(klon_glo) real :: oceanarea_glo real :: cell_area_glo(klon_glo) !#else ! real :: runoffamount_glo(ngrid) ! real :: zqsurf_iliq_glo(ngrid) !#endif if(firstcall)then oceanbulkavg=.false. oceanalbvary=.false. write(*,*)"Activate runnoff into oceans?" activerunoff=.false. call getin_p("activerunoff",activerunoff) write(*,*)" activerunoff = ",activerunoff if (activerunoff) then ALLOCATE(runoff(ngrid)) runoff(1:ngrid)=0 endif ivap=igcm_h2o_vap iliq=igcm_h2o_vap iice=igcm_h2o_ice write(*,*) "hydrol: ivap=",ivap write(*,*) " iliq=",iliq write(*,*) " iice=",iice ! Here's the deal: iice is used in place of igcm_h2o_ice both on the ! surface and in the atmosphere. ivap is used in ! place of igcm_h2o_vap ONLY in the atmosphere, while ! iliq is used in place of igcm_h2o_vap ONLY on the ! surface. ! Soon to be extended to the entire water cycle... ! LOCAL ocean surface area oceanarea=0. do ig=1,ngrid if(nint(rnat(ig)).eq.0)then oceanarea=oceanarea+cell_area(ig) endif enddo if(oceanbulkavg.and.(oceanarea.le.0.))then print*,'How are we supposed to average the ocean' print*,'temperature, when there are no oceans?' call abort endif if(activerunoff.and.(oceanarea.le.0.))then print*,'You have enabled runoff, but you have no oceans.' print*,'Where did you think the water was going to go?' call abort endif firstcall = .false. endif ! write (*,*) "oceanarea", oceanarea ! add physical tendencies already calculated ! ------------------------------------------ do ig=1,ngrid ztsurf(ig) = tsurf(ig) + ptimestep*pdtsurf(ig) pdtsurf_hyd(ig)=0.0 do iq=1,nq zqsurf(ig,iq) = qsurf(ig,iq) + ptimestep*dqsurf(ig,iq) enddo enddo do ig=1,ngrid do iq=1,nq dqs_hyd(ig,iq) = 0.0 enddo enddo do ig = 1, ngrid ! Ocean regions (rnat = 0) ! ----------------------- if(nint(rnat(ig)).eq.0)then ! Parameterization (not used for the moment) to compute the effect of solar zenith angle on the albedo ! -------------------------- ! ! if(diurnal.and.oceanalbvary)then ! fauxo = ( 1.47 - ACOS( mu0(ig) ) )/0.15 ! where does this come from (Benjamin)? ! albedo(ig) = 1.1*( .03 + .630/( 1. + fauxo*fauxo)) ! albedo(ig) = MAX(MIN(albedo(ig),0.60),0.04) ! else ! ! do nw=1,L_NSPECTV ! albedo(ig,nw) = alb_ocean ! For now, alb_ocean is defined in inifis_mod.F90. Later we could introduce spectral dependency for alb_ocean. ! enddo ! end if ! we first start by fixing the albedo of oceanic grid to that of the ocean do nw=1,L_NSPECTV albedo(ig,nw) = alb_ocean ! For now, alb_ocean is defined in inifis_mod.F90. Later we could introduce spectral dependency for alb_ocean. enddo if(ok_slab_ocean) then ! if ocean heat transport param activated frac_snow = MAX(0.0,MIN(1.0,zqsurf(ig,iice)/snow_min)) ! Critical snow height (in kg/m2) from ocean_slab_ice routine. ! Standard value should be 15kg/m2 (i.e. about 5 cm). Note that in the previous ocean param. (from BC2014), this value was 45kg/m2 (i.e. about 15cm). ! Albedo final calculation : do nw=1,L_NSPECTV alb_ice=albedo_snow_SPECTV(nw)-(albedo_snow_SPECTV(nw)-alb_ice_min)*exp(-sea_ice(ig)/h_alb_ice) ! this replaces the formulation from BC2014 ! More details on the parameterization of sea ice albedo vs thickness is provided in the wiki : ! https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Slab_ocean_model ! sea_ice is the ice thickness (calculated in ocean_slab routine) in kg/m2 ; h_alb_ice is fixed to 275.1kg/m2 i.e. 30cm based on comparisons with Brandt et al. 2005 albedo(ig,nw) = pctsrf_sic(ig)* & (albedo_snow_SPECTV(nw)*frac_snow + alb_ice*(1.0-frac_snow)) & + (1.-pctsrf_sic(ig))*alb_ocean enddo ! Oceanic ice height, just for diagnostics hice(ig) = MIN(10.,sea_ice(ig)/rhowater) else !ok_slab_ocean ; here this is the case where we are dealing with a static ocean ! calculate oceanic ice height including the latent heat of ice formation ! hice is the height of oceanic ice with a maximum of maxicethick. hice(ig) = zqsurf(ig,iice)/rhowater ! update hice to include recent snowfall twater(ig) = ztsurf(ig) - hice(ig)*RLFTT*rhowater/pcapcal(ig) ! this is the temperature water would have if we melted the entire ocean ice layer hicebis(ig) = hice(ig) hice(ig) = 0. if(twater(ig) .lt. T_h2O_ice_liq)then E=min((T_h2O_ice_liq+Tsaldiff-twater(ig))*pcapcal(ig),RLFTT*rhowater*maxicethick) hice(ig) = E/(RLFTT*rhowater) hice(ig) = max(hice(ig),0.0) hice(ig) = min(hice(ig),maxicethick) pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))*RLFTT*rhowater/pcapcal(ig)/ptimestep do nw=1,L_NSPECTV albedo(ig,nw) = albedo_snow_SPECTV(nw) ! Albedo of ice has been replaced by albedo of snow here. MT2015. enddo ! if (zqsurf(ig,iice).ge.snowlayer) then ! albedo(ig) = albedoice ! else ! albedo(ig) = albedoocean & ! + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer ! endif else pdtsurf_hyd(ig) = -hicebis(ig)*RLFTT*rhowater/pcapcal(ig)/ptimestep DO nw=1,L_NSPECTV albedo(ig,nw) = alb_ocean ENDDO endif zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice)) zqsurf(ig,iice) = hice(ig)*rhowater endif!(ok_slab_ocean) ! Continental regions (rnat = 1) ! ----------------------- elseif (nint(rnat(ig)).eq.1) then ! melt the snow if(ztsurf(ig).gt.T_h2O_ice_liq)then if(zqsurf(ig,iice).gt.1.0e-8)then a = (ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT b = zqsurf(ig,iice) fsnoi = min(a,b) zqsurf(ig,iice) = zqsurf(ig,iice) - fsnoi zqsurf(ig,iliq) = zqsurf(ig,iliq) + fsnoi ! thermal effects pdtsurf_hyd(ig) = -fsnoi*RLFTT/pcapcal(ig)/ptimestep endif else ! freeze the water if(zqsurf(ig,iliq).gt.1.0e-8)then a = -(ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT b = zqsurf(ig,iliq) fsnoi = min(a,b) zqsurf(ig,iice) = zqsurf(ig,iice) + fsnoi zqsurf(ig,iliq) = zqsurf(ig,iliq) - fsnoi ! thermal effects pdtsurf_hyd(ig) = +fsnoi*RLFTT/pcapcal(ig)/ptimestep endif endif ! add runoff (to simulate transport of water from continental regions to oceanic regions ; in practice, this prevents liquid water to build up too much on continental regions) if(activerunoff)then runoff(ig) = max(zqsurf(ig,iliq) - mx_eau_sol, 0.0) if(ngrid.gt.1)then ! runoff only exists in 3D if(runoff(ig).ne.0.0)then zqsurf(ig,iliq) = mx_eau_sol ! note: runoff is added to ocean at end endif end if endif ! re-calculate continental albedo DO nw=1,L_NSPECTV albedo(ig,nw) = albedo_bareground(ig) ENDDO if (zqsurf(ig,iice).ge.snowlayer) then DO nw=1,L_NSPECTV albedo(ig,nw) = albedo_snow_SPECTV(nw) ENDDO else DO nw=1,L_NSPECTV albedo(ig,nw) = albedo_bareground(ig) & + (albedo_snow_SPECTV(nw) - albedo_bareground(ig)) & *zqsurf(ig,iice)/snowlayer ENDDO endif else print*,'Surface type not recognised in hydrol.F!' print*,'Exiting...' call abort endif end do ! ig=1,ngrid ! simple parameterization to perform crude bulk averaging of temperature in ocean ! ---------------------------------------------------- if(oceanbulkavg)then oceanarea2=0. DO ig=1,ngrid if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then oceanarea2=oceanarea2+cell_area(ig)*pcapcal(ig) end if END DO tsea=0. DO ig=1,ngrid if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then tsea=tsea+ztsurf(ig)*cell_area(ig)*pcapcal(ig)/oceanarea2 end if END DO DO ig=1,ngrid if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0))then pdtsurf_hyd(ig) = pdtsurf_hyd(ig) + (tsea-ztsurf(ig))/oceantime end if END DO print*,'Mean ocean temperature = ',tsea,' K' endif ! shove all the runoff water into the ocean ! ----------------------------------------- if(activerunoff)then ! totalrunoff=0. do ig=1,ngrid runoffamount(ig) = cell_area(ig)*runoff(ig) ! if (nint(rnat(ig)).eq.1) then ! totalrunoff = totalrunoff + cell_area(ig)*runoff(ig) ! endif enddo ! collect on the full grid call gather(runoffamount,runoffamount_glo) call gather(zqsurf(1:ngrid,iliq),zqsurf_iliq_glo) call gather(rnat,rnat_glo) call gather(cell_area,cell_area_glo) if (is_master) then totalrunoff=0. oceanarea_glo=0. do ig=1,klon_glo if (nint(rnat_glo(ig)).eq.1) then totalrunoff = totalrunoff + runoffamount_glo(ig) endif if (nint(rnat_glo(ig)).eq.0) then oceanarea_glo = oceanarea_glo + cell_area_glo(ig) endif enddo do ig=1,klon_glo if (nint(rnat_glo(ig)).eq.0) then zqsurf_iliq_glo(ig) = zqsurf_iliq_glo(ig) + & totalrunoff/oceanarea_glo endif enddo endif! is_master ! scatter the field back on all processes call scatter(zqsurf_iliq_glo,zqsurf(1:ngrid,iliq)) ! do ig=1,ngrid ! if (nint(rnat(ig)).eq.0) then ! zqsurf(ig,iliq) = zqsurf(ig,iliq) + & ! totalrunoff/oceanarea ! endif ! enddo endif !activerunoff ! Re-add the albedo effects of CO2 ice if necessary ! ------------------------------------------------- if(co2cond)then do ig=1,ngrid if (qsurf(ig,igcm_co2_ice).gt.1.) then ! Condition changed - Need now ~1 mm CO2 ice coverage. MT2015 DO nw=1,L_NSPECTV albedo(ig,nw) = albedo_co2_ice_SPECTV(nw) ENDDO endif enddo ! ngrid endif ! co2cond do ig=1,ngrid ! We calculate here the tracer tendencies. Don't forget that we have to retrieve the dqsurf tendencies we added at the beginning of the routine ! dqs_hyd(ig,iliq)=(zqsurf(ig,iliq) - qsurf(ig,iliq))/ptimestep - dqsurf(ig,iliq) dqs_hyd(ig,iice)=(zqsurf(ig,iice) - qsurf(ig,iice))/ptimestep - dqsurf(ig,iice) enddo if (activerunoff) then call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff) endif return end subroutine hydrol