subroutine surfprop(ngrid,nq,pfra,qsurf,tsurface, & tidat,capcal,adjust,dist,albedo,emis,fluold, & ptimestep,zls) use comgeomfi_h use comsoil_h implicit none !================================================================== ! Purpose ! ------- ! set up the properties of pluto's surface ! ! Inputs ! ------ ! ngrid Number of vertical columns ! nlayer Number of layers ! qsurf(ngrid,iq) Amount of ice on surface kg/m2 ! tsurface(ngrid) surface temperature ! ! Outputs ! ------- ! albedo(ngrid) ! emis(ngrid) ! ! Both ! ---- ! ! Authors ! ------- ! Tanguy Bertrand ! !================================================================== #include "dimensions.h" #include "dimphys.h" #include "surfdat.h" #include "callkeys.h" #include "tracer.h" #include "comcstfi.h" !----------------------------------------------------------------------- ! Arguments INTEGER ngrid, nq REAL,INTENT(IN) :: qsurf(ngrid,nq) REAL,INTENT(IN) :: fluold(ngrid,nq) REAL,INTENT(IN) :: ptimestep REAL,INTENT(IN) :: zls REAL,INTENT(IN) :: tsurface(ngrid) REAL,INTENT(IN) :: capcal(ngrid) REAL,INTENT(IN) :: adjust REAL,INTENT(IN) :: dist REAL,INTENT(OUT) :: albedo(ngrid) REAL,INTENT(OUT) :: emis(ngrid) REAL,INTENT(OUT) :: tidat(ngrid,nsoilmx) REAL,INTENT(IN) :: pfra(ngrid) !----------------------------------------------------------------------- ! Local variables REAL stephan DATA stephan/5.67e-08/ ! Stephan Boltzman constant SAVE stephan REAL tsa,tsb,coef,aa,bb REAL emin,emax,emax2,tid REAL n2cover,n2coversun,gamm,frac,pls,facls INTEGER ig,isoil,ice_case LOGICAL firstcall SAVE firstcall DATA firstcall/.true./ !----------------------------------------------------------------------- ! Special aging for methane (not implemented) REAL albstep ! Time constant for albedo change (age) DATA albstep/1.e7/ REAL albmet(ngridmx) SAVE albmet REAL albmetmin ! Min alb for methane ice DATA albmetmin/0.5/ REAL albmetmax ! Max alb for methane ice DATA albmetmax/0.8/ !----------------------------------------------------------------------- ! 1) ALBEDOS and EMISSIVITY ! A. N2 ! CASE (0) ! fixed albedo ! CASE (1) ! Albedo decreases with thickness ! CASE (2) ! Special Sputnik differences of albedo ! CASE (3) ! Albedo increases (delta neg) or decreases (delta pos) with sublimationi rates ! CASE (4) ! Albedo Difference in N/S (e.g. used for Triton) ! CASE (5) ! Special Sputnik differences of albedo in small (1 pixel) patches (e.g. simulating dark patches / plumes) ! --> EMISSIVITY N2: based on the alpha/beta transition ! B. CO ! C. CH4 ! CASE (0) ! 2 albedos, one for the tropics, one for the poles ! CASE (1) ! 3 albedos, one for the tropics, 2 for the poles (north and south) ! CASE (2) ! 2 albedos + albedo feedback ! SELECT CASE (feedback_met) ! CASE (0) ! Default (linear from alb_ch4_eq) ! CASE (1) ! Hyperbolic tangent old ! CASE (2) ! hyperbolic tangent old ! CASE (3) ! hyperbolic tangent equation with parameters ! CASE (3) ! Eq, poles N, pole S + depending on Ls ! D. Tholins ! CASE (0) ! Default, 2 albedos, one for the tropics, one for the poles ! CASE (1) ! Special mode one region with a different albedo ! E. Tholins read from file ! specalb ! 2) Changes of Thermal inertia with time ! if (changeti) : change of seasonal TI ! if (changetid) : change of diurnal TI ! 3) OTHER TESTS !----------------------------------------------------------------------- ! 1) ALBEDOS and EMISSIVITY !----------------------------------------------------------------------- ! Some input parameters: pls=zls*180./pi ! for equation of feedbackalbedo if (methane.and.mode_ch4.eq.2.and.feedback_met.eq.3) then aa=fdch4_finalb-fdch4_depalb bb=fdch4_finalb endif ! Loop on all points DO ig=1,ngrid ! Looking for dominant ice: ice_case=0 ! Tholins if (methane) then if (qsurf(ig,igcm_ch4_ice).ge.thres_ch4ice) then ice_case=1 endif endif if (carbox) then if (qsurf(ig,igcm_co_ice).ge.thres_coice) then ice_case=2 endif endif if (qsurf(ig,igcm_n2).ge.thres_n2ice) then ice_case=3 endif !--------------------------------- ! 1.A. N2 !--------------------------------- if (ice_case.eq.3) then ! EMISSIVITY N2 ! change emis if phase alpha different de phase beta frac=0.5*(1.+tanh(6.*(35.6-tsurface(ig))/0.5)) if (frac.gt.1.) frac=1. emis(ig)=frac*emis_n2a+(1.-frac)*emis_n2b ! ALBEDO N2 and emissivity changes SELECT CASE (mode_n2) CASE (0) ! fixed albedo albedo(ig)=min(max(alb_n2b+adjust,0.),1.) CASE (1) ! Albedo decreases with thickness albedo(ig)=(alb_n2b-deltab)/(1.-10000.)*qsurf(ig,igcm_n2)+alb_n2b albedo(ig)=min(max(alb_n2b-deltab,albedo(ig)),alb_n2b) ! should not be higher than alb_n2b and not lower than alb_n2b-deltab CASE (2) ! Special Sputnik differences of albedo if (qsurf(ig,igcm_n2).ge.1.e4.and.(long(ig)*180./pi.le.-161.or.long(ig)*180./pi.ge.146)) then if ( (lati(ig)*180./pi.ge.25.).or. & ((long(ig)*180./pi.gt.140.).and. & (long(ig)*180./pi.lt.165.)) ) then albedo(ig)=alb_n2b-deltab else albedo(ig)=alb_n2b+deltab endif else albedo(ig)=alb_n2b endif CASE (3) ! Albedo increases (delta neg) or decreases (delta pos) with sublimation rates albedo(ig)=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig) albedo(ig)=min(max(alb_n2b-deltab,albedo(ig)),alb_n2b+deltab) !! Triton: !albedo(ig)=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig) !albedo(ig)=min(max(alb_n2b-deltab,albedo(ig)),alb_n2b+deltab) CASE (4) ! Albedo Difference in N/S if (lati(ig)*180./pi.ge.0.) then albedo(ig)=min(max(alb_n2b-deltab+adjust,0.),1.) else albedo(ig)=min(max(alb_n2b+deltab+adjust,0.),1.) endif CASE (5) ! Special Sputnik differences of albedo in patches !! Patches Nogcm !if (qsurf(ig,igcm_n2).ge.1.e4.and.(long(ig)*180./pi.le.180.).and.(long(ig)*180./pi.ge.174.) ) then ! if (((lati(ig)*180./pi.le.46.).and.(lati(ig)*180./pi.ge.42.)) & ! .or. ((lati(ig)*180./pi.le.36.).and.(lati(ig)*180./pi.ge.32.)) & ! .or. ((lati(ig)*180./pi.le.26.).and.(lati(ig)*180./pi.ge.22.)) & ! .or. ((lati(ig)*180./pi.le.16.).and.(lati(ig)*180./pi.ge.12.)) & ! .or. ((lati(ig)*180./pi.le.6.).and.(lati(ig)*180./pi.ge.2.)) & ! ) then ! albedo(ig)=alb_n2b-deltab !! Patches GCM if (qsurf(ig,igcm_n2).ge.1.e4) then if (((lati(ig)*180./pi.lt.33.).and.(lati(ig)*180./pi.gt.32.).and. & (long(ig)*180./pi.lt.165.5).and.(long(ig)*180./pi.gt.164.5)) & .or. ((lati(ig)*180./pi.lt.30.5).and.(lati(ig)*180./pi.gt.29.5).and. & (long(ig)*180./pi.lt.169.).and.(long(ig)*180./pi.gt.168.))) then albedo(ig)=alb_n2b-deltab else if (((lati(ig)*180./pi.lt.30.5).and.(lati(ig)*180./pi.gt.29.).and. & (long(ig)*180./pi.lt.165.5).and.(long(ig)*180./pi.gt.164.5)) & .or. ((lati(ig)*180./pi.lt.33.).and.(lati(ig)*180./pi.gt.32.).and. & (long(ig)*180./pi.lt.169.).and.(long(ig)*180./pi.gt.168.))) then albedo(ig)=alb_n2b+deltab else albedo(ig)=alb_n2b endif else albedo(ig)=alb_n2b endif CASE (7) ! Albedo from albedodat and adjusted emissivity albedo(ig)=albedodat(ig) ! adjust emissivity if convergeps = true emis(ig)=min(max(emis(ig)+adjust,0.),1.) CASE DEFAULT write(*,*) 'STOP in surfprop.F90: mod_n2 not found' stop END SELECT !--------------------------------- ! 1.B. CO !--------------------------------- else if (ice_case.eq.2) then albedo(ig)=alb_co emis(ig)=emis_co !--------------------------------- ! 1.C. CH4 !--------------------------------- else if (ice_case.eq.1) then SELECT CASE (mode_ch4) CASE (0) ! 2 albedos, one for the tropics, one for the poles emis(ig)=emis_ch4 if (lati(ig)*180./pi.le.metlateq.and.lati(ig)*180./pi.gt.-metlateq) then albedo(ig)=alb_ch4_eq else albedo(ig)=alb_ch4 endif CASE (1) ! 3 albedos, one for the tropics, 2 for the poles (north and south) emis(ig)=emis_ch4 if (lati(ig)*180./pi.le.metlateq.and.lati(ig)*180./pi.gt.-metlateq) then albedo(ig)=alb_ch4_eq else if (lati(ig)*180./pi.le.-metlateq) then albedo(ig)=alb_ch4_s else albedo(ig)=alb_ch4 endif CASE (2) ! 2 albedos + albedo feedback emis(ig)=emis_ch4 albedo(ig)=alb_ch4 if (lati(ig)*180./pi.le.metlateq.and.lati(ig)*180./pi.gt.-metlateq) then albedo(ig)=alb_ch4_eq endif !! Albedo feedback if (lati(ig)*180./pi.le.fdch4_latn.and.lati(ig)*180./pi.gt.fdch4_lats) then if (long(ig)*180./pi.le.fdch4_lone.and.long(ig)*180./pi.gt.fdch4_lonw) then if (qsurf(ig,igcm_ch4_ice).lt.fdch4_maxice) then ! fd would not apply on BTD SELECT CASE (feedback_met) CASE (0) ! Default (linear from alb_ch4_eq) albedo(ig)=(1.-alb_ch4_eq)/0.002*qsurf(ig,igcm_ch4_ice)+alb_ch4_eq !emis(ig)=(1.-emis_ch4)/0.002*qsurf(ig,igcm_ch4_ice)+emis_ch4 if (albedo(ig).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb !if (emis(ig).gt.1.) emis(ig)=1. CASE (1) ! Hyperbolic tangent old albedo(ig)=0.6*0.5*(1.+tanh(6.*(0.001+qsurf(ig,igcm_ch4_ice))/0.5))+0.3 ! !emis(ig)=0.2*0.5*(1.+tanh(6.*(0.001+qsurf(ig,igcm_ch4_ice))/0.5))+0.8 ! if (albedo(ig).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb !if (emis(ig).gt.1.) emis(ig)=1. CASE (2) ! hyperbolic tangent old albedo(ig)=0.5*0.6*(1.+tanh(1000.*(qsurf(ig,igcm_ch4_ice)-0.002)))+0.3 ! !emis(ig)=0.2*0.5*(1.+tanh(1000.*(qsurf(ig,igcm_ch4_ice)-0.002)))+0.8 ! if (albedo(ig).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb !if (emis(ig).gt.1.) emis(ig)=1. CASE (3) ! hyperbolic tangent equation with parameters albedo(ig)=aa*(-1+tanh(fdch4_ampl*qsurf(ig,igcm_ch4_ice)))+bb if (albedo(ig).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb CASE DEFAULT write(*,*) 'STOP surfprop.F90: fd wrong' stop END SELECT endif endif endif CASE (3) ! Eq, poles N, pole S + depending on Ls emis(ig)=emis_ch4 if (lati(ig)*180./pi.le.metlateq.and.lati(ig)*180./pi.gt.-metlateq) then albedo(ig)=alb_ch4_eq ! Pure methane ice else if (lati(ig)*180./pi.le.-metlateq) then albedo(ig)=alb_ch4_s ! Pure methane ice if (pls.le.metls2.and.pls.gt.metls1) then albedo(ig)=alb_ch4_s+deltab ! Also used for N2, careful endif else albedo(ig)=alb_ch4 ! Pure methane ice endif CASE (4) ! Albedo from albedodat emis(ig)=emis_ch4 albedo(ig)=albedodat(ig) CASE DEFAULT write(*,*) 'STOP in surfprop.F90:mod_ch4 not found' stop END SELECT !--------------------------------- ! 1.D. THOLINS !--------------------------------- else SELECT CASE (mode_tholins) CASE (0) ! Default, 2 albedos, one for the tropics, one for the poles if (lati(ig)*180./pi.le.tholateq.and.lati(ig)*180./pi.gt.-tholateq) then albedo(ig)=alb_tho_eq emis(ig)=emis_tho_eq else albedo(ig)=alb_tho ! default = 0.1 emis(ig)=emis_tho endif CASE (1) ! Special mode one region with a different albedo if (lati(ig)*180./pi.le.tholateq.and.lati(ig)*180./pi.gt.-tholateq) then albedo(ig)=alb_tho_eq emis(ig)=emis_tho_eq else albedo(ig)=alb_tho ! default = 0.1 emis(ig)=emis_tho endif if (lati(ig)*180./pi.le.tholatn.and.lati(ig)*180./pi.ge.tholats & .and.long(ig)*180./pi.ge.tholone.and.long(ig)*180./pi.le.tholonw) then albedo(ig)=alb_tho_spe emis(ig)=emis_tho_spe endif CASE (2) ! Depends on the albedo map albedo(ig)=albedodat(ig) if (albedo(ig).gt.alb_ch4) then emis(ig)=emis_ch4 else emis(ig)=emis_tho endif CASE DEFAULT write(*,*) 'STOP in surfprop.F90:mod_ch4 not found' stop END SELECT !--------------------------------- ! 1.E. Tholins read from file !--------------------------------- if (specalb) then albedo(ig)=albedodat(ig) ! Specific ground properties !if (albedodat(ig).lt.0.25) then ! albedo(ig)=alb_tho_eq !else if (albedodat(ig).lt.0.40) then ! albedo(ig)=0.35 !else if (albedodat(ig).lt.0.70) then ! albedo(ig)=0.51 !endif endif endif ! ice_case ENDDO ! ig ngrid !----------------------------------------------------------------------- ! 2) Changes of Thermal inertia with time !----------------------------------------------------------------------- IF (changeti) then ! change of seasonal TI facls=8. DO ig=1,ngrid ! get depth of the different layers ! limit diurnal / seasonal if (changetid) then if (qsurf(ig,igcm_n2)>1.e-3) then emin=facls*ITN2d/volcapa*(daysec/pi)**0.5 tid=ITN2d else if (qsurf(ig,igcm_ch4_ice)>1.e-3) then emin=facls*ITCH4d/volcapa*(daysec/pi)**0.5 tid=ITCH4d else emin=facls*ITH2Od/volcapa*(daysec/pi)**0.5 tid=ITH2Od endif else if (ngrid.ne.1) then emin=facls*inertiedat(ig,1)/volcapa*(daysec/pi)**0.5 tid=inertiedat(ig,1) else emin=facls*ITH2Od/volcapa*(daysec/pi)**0.5 tid=ITH2Od endif ! limit for N2 emax=emin+qsurf(ig,igcm_n2)/1000. ! limit for CH4 if (methane) then emax2=emax+qsurf(ig,igcm_ch4_ice)/1000. else emax2=0. endif do isoil=0,nsoilmx-1 if (mlayer(isoil).le.emin) then ! diurnal TI tidat(ig,isoil+1)=tid else if (isoil.gt.0.and.(mlayer(isoil).gt.emin).and.(mlayer(isoil-1).lt.emin)) then ! still diurnal TI tidat(ig,isoil+1)=tid else if ((mlayer(isoil).gt.emin).and.(mlayer(isoil).le.emax)) then ! TI N2 tidat(ig,isoil+1)=ITN2 else if ((mlayer(isoil).gt.emin).and.(mlayer(isoil).le.emax2)) then tidat(ig,isoil+1)=ITCH4 else tidat(ig,isoil+1)=ITH2O endif enddo ENDDO print*, emin print*, tidat(1,:) print*, mlayer(:) ELSE DO ig=1,ngrid tidat(ig,:)=inertiedat(ig,:) ENDDO ENDIF !----------------------------------------------------------------------- ! 3) Tests changements emis !----------------------------------------------------------------------- !n2cover=0. !n2coversun=0. !DO ig=1,ngrid ! if (qsurf(ig,igcm_n2).ge.0.001) then ! n2cover=n2cover+area(ig) ! if (pfra(ig).gt.0.) then ! n2coversun=n2coversun+area(ig) ! endif ! endif !enddo !gamm=n2cover/n2coversun !gamm=1. !tsb=(1/gamm*Fat1AU/dist**2*(1.-alb_n2b)/(emis_n2b*stephan))**(1/4.) !tsa=(1/gamm*Fat1AU/dist**2*(1.-alb_n2b)/(emis_n2a*stephan))**(1/4.) !frac=max((35.6-tsb)/abs(tsa-tsb),0.) !write(*,*) 'gamm=',gamm,n2cover,n2coversun !write(*,*) 'tsb=',tsb !write(*,*) 'tsa=',tsa !write(55,*) 'frac=',frac,tsb,tsa end subroutine surfprop