Changeset 3329
- Timestamp:
- May 15, 2024, 8:29:23 PM (7 months ago)
- Location:
- trunk
- Files:
-
- 6 added
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.PLUTO.old/libf/phypluto/physiq.F
r3275 r3329 701 701 !global mean 1D radiative equiilium 702 702 ig=1 703 mu0(ig) = 2**(-0.5)703 mu0(ig) = 0.5 704 704 fract(ig)= 1/(4*mu0(ig)) 705 705 !print*,'WARNING GLOBAL MEAN CALCULATION' -
trunk/LMDZ.PLUTO/deftank/callphys.def
r3197 r3329 88 88 89 89 naerkind = 1 90 hazeprop = optprop_tholins_fractal100nm.dat90 hazeprop_file = optprop_tholins_fractal100nm.dat 91 91 92 92 ## Other physics options -
trunk/LMDZ.PLUTO/deftank/z2sig.def
r3193 r3329 1 1 18.00000 2 0.010 3 0.015 2 0.010 3 0.015 4 4 0.025 5 5 0.040 -
trunk/LMDZ.PLUTO/libf/phypluto/aerosol_mod.F90
r3195 r3329 74 74 IF (firstcall) then 75 75 firstcall=.false. 76 file_path=trim(datadir)//'/haze_prop/ hazemmr.txt'76 file_path=trim(datadir)//'/haze_prop/'//hazemmr_file 77 77 open(224,file=file_path,form='formatted') 78 78 do ifine=1,Nfine -
trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90
r3275 r3329 7 7 subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, & 8 8 albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & 9 zzlay,tsurf,fract,dist_star,aerosol,muvar, 9 zzlay,tsurf,fract,dist_star,aerosol,muvar, & 10 10 dtlw,dtsw,fluxsurf_lw, & 11 11 fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw, & … … 26 26 use ioipsl_getin_p_mod, only: getin_p 27 27 use gases_h, only: ngasmx 28 use radii_mod, only : su_aer_radii 28 use radii_mod, only : su_aer_radii, haze_reffrad_fix 29 29 use aerosol_mod, only : iaero_haze 30 30 use aeropacity_mod, only: aeropacity … … 225 225 parameter(Nfine=701) 226 226 real,save :: levdat(Nfine),vmrdat(Nfine) 227 REAL dtlw_hcn_c2h2(ngrid, nlayer) ! cooling rate (K/s) due to C2H2/HCN (diagnostic) 227 228 real :: vmrch4(ngrid,nlayer) ! vmr ch4 from vmrch4_proffix 228 229 … … 506 507 507 508 if (aerohaze) then 508 !if (haze_radproffix) then509 !call haze_reffrad_fix(ngrid,nlayer,zzlay,&510 !reffrad,nueffrad)511 !endif509 if (haze_radproffix) then 510 call haze_reffrad_fix(ngrid,nlayer,zzlay,& 511 reffrad,nueffrad) 512 endif 512 513 513 514 ! Get 3D aerosol optical properties. … … 554 555 ! dtlw_hcn_c2h2=0. 555 556 if (cooling) then 556 !CALL cooling_hcn_c2h2(ngrid,nlayer,pplay,&557 !pt,dtlw_hcn_c2h2)557 CALL cooling_hcn_c2h2(ngrid,nlayer,pplay,& 558 pt,dtlw_hcn_c2h2) 558 559 endif 559 560 … … 575 576 576 577 577 if (aerohaze) then 578 do iaer=1,naerkind 579 ! Shortwave. 580 do nw=1,L_NSPECTV 581 582 do l=1,nlayer 583 584 temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer) & 585 *QREFvis3d(ig,nlayer+1-l,iaer) 586 587 temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer) & 588 *QREFvis3d(ig,max(nlayer-l,1),iaer) 589 590 qxvaer(2*l,nw,iaer) = temp1 591 qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 592 593 temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer) 594 temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer) 595 596 qsvaer(2*l,nw,iaer) = temp1 597 qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 598 599 temp1=gvis3d(ig,nlayer+1-l,nw,iaer) 600 temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer) 601 602 gvaer(2*l,nw,iaer) = temp1 603 gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 604 605 end do ! nlayer 606 607 qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer) 608 qxvaer(2*nlayer+1,nw,iaer)=0. 609 610 qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer) 611 qsvaer(2*nlayer+1,nw,iaer)=0. 612 613 gvaer(1,nw,iaer)=gvaer(2,nw,iaer) 614 gvaer(2*nlayer+1,nw,iaer)=0. 615 616 end do ! L_NSPECTV 617 618 do nw=1,L_NSPECTI 619 ! Longwave 620 do l=1,nlayer 621 622 temp1=QIRsQREF3d(ig,nlayer+1-l,nw,iaer) & 623 *QREFir3d(ig,nlayer+1-l,iaer) 624 625 temp2=QIRsQREF3d(ig,max(nlayer-l,1),nw,iaer) & 626 *QREFir3d(ig,max(nlayer-l,1),iaer) 627 628 qxiaer(2*l,nw,iaer) = temp1 629 qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2 630 631 temp1=temp1*omegair3d(ig,nlayer+1-l,nw,iaer) 632 temp2=temp2*omegair3d(ig,max(nlayer-l,1),nw,iaer) 633 634 qsiaer(2*l,nw,iaer) = temp1 635 qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2 636 637 temp1=gir3d(ig,nlayer+1-l,nw,iaer) 638 temp2=gir3d(ig,max(nlayer-l,1),nw,iaer) 639 640 giaer(2*l,nw,iaer) = temp1 641 giaer(2*l+1,nw,iaer)=(temp1+temp2)/2 642 643 end do ! nlayer 644 645 qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer) 646 qxiaer(2*nlayer+1,nw,iaer)=0. 647 648 qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer) 649 qsiaer(2*nlayer+1,nw,iaer)=0. 650 651 giaer(1,nw,iaer)=giaer(2,nw,iaer) 652 giaer(2*nlayer+1,nw,iaer)=0. 653 654 end do ! L_NSPECTI 655 656 end do ! naerkind 657 658 ! Test / Correct for freaky s. s. albedo values. 659 do iaer=1,naerkind 660 do k=1,L_LEVELS 661 662 do nw=1,L_NSPECTV 663 if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then 664 message='Serious problems with qsvaer values' 665 call abort_physic(subname,message,1) 666 endif 667 if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then 668 qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer) 669 endif 670 end do 671 672 do nw=1,L_NSPECTI 673 if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then 674 message='Serious problems with qsvaer values' 675 call abort_physic(subname,message,1) 676 endif 677 if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then 678 qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer) 679 endif 680 end do 681 682 end do ! L_LEVELS 683 end do ! naerkind 684 685 endif ! aerohaze 578 do iaer=1,naerkind 579 ! Shortwave. 580 do nw=1,L_NSPECTV 581 582 do l=1,nlayer 583 584 temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer) & 585 *QREFvis3d(ig,nlayer+1-l,iaer) 586 587 temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer) & 588 *QREFvis3d(ig,max(nlayer-l,1),iaer) 589 590 qxvaer(2*l,nw,iaer) = temp1 591 qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 592 593 temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer) 594 temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer) 595 596 qsvaer(2*l,nw,iaer) = temp1 597 qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 598 599 temp1=gvis3d(ig,nlayer+1-l,nw,iaer) 600 temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer) 601 602 gvaer(2*l,nw,iaer) = temp1 603 gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 604 605 end do ! nlayer 606 607 qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer) 608 qxvaer(2*nlayer+1,nw,iaer)=0. 609 610 qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer) 611 qsvaer(2*nlayer+1,nw,iaer)=0. 612 613 gvaer(1,nw,iaer)=gvaer(2,nw,iaer) 614 gvaer(2*nlayer+1,nw,iaer)=0. 615 616 end do ! L_NSPECTV 617 618 do nw=1,L_NSPECTI 619 ! Longwave 620 do l=1,nlayer 621 622 temp1=QIRsQREF3d(ig,nlayer+1-l,nw,iaer) & 623 *QREFir3d(ig,nlayer+1-l,iaer) 624 625 temp2=QIRsQREF3d(ig,max(nlayer-l,1),nw,iaer) & 626 *QREFir3d(ig,max(nlayer-l,1),iaer) 627 628 qxiaer(2*l,nw,iaer) = temp1 629 qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2 630 631 temp1=temp1*omegair3d(ig,nlayer+1-l,nw,iaer) 632 temp2=temp2*omegair3d(ig,max(nlayer-l,1),nw,iaer) 633 634 qsiaer(2*l,nw,iaer) = temp1 635 qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2 636 637 temp1=gir3d(ig,nlayer+1-l,nw,iaer) 638 temp2=gir3d(ig,max(nlayer-l,1),nw,iaer) 639 640 giaer(2*l,nw,iaer) = temp1 641 giaer(2*l+1,nw,iaer)=(temp1+temp2)/2 642 643 end do ! nlayer 644 645 qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer) 646 qxiaer(2*nlayer+1,nw,iaer)=0. 647 648 qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer) 649 qsiaer(2*nlayer+1,nw,iaer)=0. 650 651 giaer(1,nw,iaer)=giaer(2,nw,iaer) 652 giaer(2*nlayer+1,nw,iaer)=0. 653 654 end do ! L_NSPECTI 655 656 end do ! naerkind 657 658 ! Test / Correct for freaky s. s. albedo values. 659 do iaer=1,naerkind 660 do k=1,L_LEVELS 661 662 do nw=1,L_NSPECTV 663 if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then 664 message='Serious problems with qsvaer values' 665 call abort_physic(subname,message,1) 666 endif 667 if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then 668 qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer) 669 endif 670 end do 671 672 do nw=1,L_NSPECTI 673 if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then 674 message='Serious problems with qsvaer values' 675 call abort_physic(subname,message,1) 676 endif 677 if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then 678 qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer) 679 endif 680 end do 681 682 end do ! L_LEVELS 683 end do ! naerkind 684 686 685 !----------------------------------------------------------------------- 687 686 ! Aerosol optical depths 688 687 !----------------------------------------------------------------------- 689 IF (aerohaze) THEN 690 do iaer=1,naerkind ! a bug was here 688 do iaer=1,naerkind ! a bug was here 691 689 do k=0,nlayer-1 692 690 … … 706 704 707 705 end do ! naerkind 708 ELSE709 tauaero(:,:)=0710 ENDIF ! aerohaze711 706 712 707 ! Albedo and Emissivity. … … 727 722 728 723 if (generic_condensation .or. methane .or. carbox) then 724 725 ! IF (methane) then 726 727 ! do l=1,nlayer 728 ! qvar(2*l) = vmrch4(ig,nlayer+1-l)/100.* & 729 ! mmol(igcm_ch4_gas)/mmol(igcm_n2) 730 ! qvar(2*l+1) = ((vmrch4(ig,nlayer+1-l)+vmrch4(ig, & 731 ! max(nlayer-l,1)))/2.)/100.* & 732 ! mmol(igcm_ch4_gas)/mmol(igcm_n2) 733 ! end do 734 ! qvar(1)=qvar(2) 735 736 ! ELSE 729 737 730 738 ! For now, only one GCS tracer can be both variable and radiatively active -
trunk/LMDZ.PLUTO/libf/phypluto/datafile_mod.F90
r3193 r3329 12 12 character(len=300),save :: datadir='/u/lmdz/WWW/planets/LMDZ.GENERIC/datagcm' 13 13 !$OMP THREADPRIVATE(datadir) 14 14 15 15 ! Subdirectories of 'datadir': 16 16 17 17 ! surfdir stores planetary topography, albedo, etc. (surface.nc files) 18 18 character(len=12),parameter :: surfdir="surface_data" 19 19 20 20 ! aerdir stores aerosol properties files (optprop_*dat files) 21 character(LEN=18),parameter :: aerdir="aerosol_properties" 21 character(LEN=18),parameter :: aerdir="aerosol_properties" 22 22 23 23 ! Data haze properties 24 character(len=300),save :: hazeprop 25 24 character(len=300),save :: hazeprop_file 25 character(len=300),save :: hazerad_file 26 character(len=300),save :: hazemmr_file 27 26 28 end module datafile_mod 27 29 !----------------------------------------------------------------------- -
trunk/LMDZ.PLUTO/libf/phypluto/dyn1d/rcm1d.F
r3258 r3329 22 22 use time_phylmdz_mod, only: daysec, dtphys, day_step, ecritphy, 23 23 & nday, iphysiq 24 use callkeys_mod, only: tracer, specOLR,pceil 24 use callkeys_mod, only: tracer, specOLR,pceil,haze 25 25 USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff, sig, 26 26 & presnivs,pseudoalt,scaleheight … … 78 78 c 79 79 INTEGER day0 ! date initial (sol ; =0 a Ls=0) 80 INTEGER lecttsoil ! lecture of tsoil from proftsoil 81 INTEGER lecthaze ! lecture of haze from profhaze 80 82 REAL day ! date durant le run 81 83 REAL time ! time (0<time<1 ; time=0.5 a midi) … … 329 331 endif 330 332 enddo !of do iq=1,nq 331 ! check for n2 / h2o_ice tracers:333 ! check for n2 / ch4_ice tracers: 332 334 i_n2=0 333 335 do iq=1,nq … … 342 344 elseif (tname(iq)=="co_gas") then 343 345 i_co_gas=iq 346 elseif (tname(iq)=="haze") then 347 i_haze=iq 344 348 endif 345 349 enddo … … 657 661 q(ilayer,iq) = 0.01*coref*28./28. 658 662 ENDDO 659 ! else if ((iq.eq.i_haze10).or.(iq.eq.i_haze30).or.(iq.eq.i_haze).or.(iq.eq.i_haze50).or.(iq.eq.i_haze100)) then 660 ! hazeref=0. ! default value for haze mmr 661 ! PRINT *,'qhaze (kg/kg) ?' 662 ! call getin("hazeref",hazeref) 663 ! write(*,*) " haze ref (kg/kg) = ",hazeref 664 ! DO ilayer=1,nlayer 665 ! q(ilayer,iq) = hazeref 666 ! ENDDO 667 ! else if (iq.eq.i_prec_haze) then 668 ! prechazeref=0. ! default value for vmr ch4 669 ! PRINT *,'qprechaze (kg/kg) ?' 670 ! call getin("prechazeref",prechazeref) 671 ! write(*,*) " prec haze ref (kg/kg) = ",prechazeref 672 ! DO ilayer=1,nlayer 673 ! q(ilayer,iq) = prechazeref 674 ! ENDDO 675 663 else if ((iq.eq.i_haze10).or.(iq.eq.i_haze30) 664 & .or.(iq.eq.i_haze).or.(iq.eq.i_haze50) 665 & .or.(iq.eq.i_haze100)) then 666 hazeref=0. ! default value for haze mmr 667 PRINT *,'qhaze (kg/kg) ?' 668 call getin("hazeref",hazeref) 669 write(*,*) " haze ref (kg/kg) = ",hazeref 670 DO ilayer=1,nlayer 671 q(ilayer,iq) = hazeref 672 ENDDO 673 else if (iq.eq.i_prec_haze) then 674 prechazeref=0. ! default value for vmr ch4 675 PRINT *,'qprechaze (kg/kg) ?' 676 call getin("prechazeref",prechazeref) 677 write(*,*) " prec haze ref (kg/kg) = ",prechazeref 678 DO ilayer=1,nlayer 679 q(ilayer,iq) = prechazeref 680 ENDDO 676 681 else 677 682 DO ilayer=1,nlayer … … 856 861 ! ------------------------------------------ 857 862 volcapa=1.e6 ! volumetric heat capacity 858 DO isoil=1,nsoil 863 lecttsoil=0 ! default value for lecttsoil 864 call getin("lecttsoil",lecttsoil) 865 if (lecttsoil==1) then 866 OPEN(14,file='proftsoil',status='old',form='formatted',err=101) 867 DO isoil=1,nsoil 868 READ (14,*) tsoil(isoil) 869 inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia 870 ENDDO 871 GOTO 201 872 101 STOP'fichier proftsoil inexistant' 873 201 CONTINUE 874 CLOSE(14) 875 876 else 877 DO isoil=1,nsoil 859 878 inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia 860 879 tsoil(isoil)=tsurf(1) ! soil temperature 861 ENDDO 880 ENDDO 881 endif 882 862 883 863 884 ! Initialize depths … … 870 891 enddo 871 892 893 ! Initialize haze profile 894 ! ------------------------------------------ 895 if (haze) then 896 897 lecthaze=0 ! default value for lecthaze 898 call getin("lecthaze",lecthaze) 899 if (lecthaze==1) then 900 OPEN(15,file='profhaze',status='old',form='formatted',err=301) 901 DO iq=1,nq 902 if (iq.eq.i_haze) then 903 DO ilayer=1,nlayer 904 READ (15,*) q(ilayer,iq) 905 ENDDO 906 endif 907 ENDDO 908 GOTO 401 909 301 STOP'fichier profhaze inexistant' 910 401 CONTINUE 911 CLOSE(15) 912 endif 913 endif 872 914 ! Initialize cloud fraction and oceanic ice !AF24: removed 873 915 ! Initialize slab ocean !AF24: removed … … 1022 1064 CLOSE(12) 1023 1065 endif 1024 1066 ! save haze profile 1067 if (haze.and.lecthaze.eq.1) then 1068 OPEN(16,file='profhaze.out',form='formatted') 1069 DO iq=1,nq 1070 if (iq.eq.i_haze) then 1071 DO ilayer=1,nlayer 1072 write(16,*) q(ilayer,iq) 1073 ENDDO 1074 endif 1075 ENDDO 1076 CLOSE(16) 1077 endif 1025 1078 1026 1079 ENDDO ! fin de la boucle temporelle -
trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90
r3275 r3329 13 13 use radcommon_h, only: ini_radcommon_h 14 14 use radii_mod, only: radfixed, Nmix_n2 15 use datafile_mod, only: datadir, hazeprop 15 use datafile_mod, only: datadir, hazeprop_file, hazerad_file, hazemmr_file 16 16 use comdiurn_h, only: sinlat, coslat, sinlon, coslon 17 17 use comgeomfi_h, only: totarea, totarea_planet … … 138 138 139 139 if (is_master) write(*,*) "Haze optical properties datafile" 140 hazeprop="../../../datadir" ! default path 141 call getin_p("hazeprop",hazeprop) 142 if (is_master) write(*,*) trim(rname)//" hazeprop = ",trim(hazeprop) 140 hazeprop_file="optprop_tholins_fractal050nm.dat" ! default file 141 call getin_p("hazeprop_file",hazeprop_file) 142 if (is_master) write(*,*) trim(rname)//" hazeprop_file = ",trim(hazeprop_file) 143 144 if (is_master) write(*,*) "Haze radii datafile" 145 hazerad_file="hazerad.txt" ! default file 146 call getin_p("hazerad_file",hazerad_file) 147 if (is_master) write(*,*) trim(rname)//" hazerad_file = ",trim(hazerad_file) 148 149 if (is_master) write(*,*) "Haze mmr datafile" 150 hazemmr_file="hazemmr.txt" ! default file 151 call getin_p("hazemmr_file",hazemmr_file) 152 if (is_master) write(*,*) trim(rname)//" hazemmr_file = ",trim(hazemmr_file) 143 153 144 154 -
trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90
r3275 r3329 54 54 fast,fasthaze,haze,metcloud,monoxcloud,& 55 55 n2cond,nearn2cond,noseason_day,conservn2, & 56 kbo,triton,paleo,paleoyears, &56 convergeps,kbo,triton,paleo,paleoyears, & 57 57 carbox, methane, oldplutovdifc, oldplutocorrk, & 58 58 aerohaze,haze_proffix,source_haze, & … … 67 67 use phys_state_var_mod 68 68 use callcorrk_mod, only: callcorrk 69 !use callcorrk_pluto_mod, only: callcorrk_pluto69 use callcorrk_pluto_mod, only: callcorrk_pluto 70 70 use vdifc_mod, only: vdifc 71 71 use vdifc_pluto_mod, only: vdifc_pluto … … 260 260 REAL,save :: tcond1p4Pa 261 261 DATA tcond1p4Pa/38/ 262 real :: tidat(ngrid,nsoilmx) ! thermal inertia soil 262 263 263 264 … … 306 307 real zlss ! Sub solar point longitude (radians). 307 308 real zday ! Date (time since Ls=0, calculated in sols). 309 REAL,save :: saveday ! saved date 310 REAL,save :: savedeclin ! saved declin 308 311 real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers. 309 312 real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries. … … 619 622 call iniorbit(apoastr,periastr,year_day,peri_day,obliquit) 620 623 624 savedeclin=0. 625 saveday=pday 626 !adjust=0. ! albedo adjustment for convergeps 621 627 622 628 ! Initialize soil. … … 871 877 if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight. 872 878 mu0 = cos(pi*szangle/180.0) 879 fract= 1/(4*mu0) ! AF24: from pluto.old 873 880 endif 874 881 875 882 endif 876 883 877 ! AF24: TODO insert surfprop for pluto & triton around here 884 885 ! Pluto albedo /IT changes depending on surface ices (only in 3D) 886 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 887 if (ngrid.ne.1) then 888 889 !! Specific to change albedo of N2 so that Psurf 890 !! converges toward 1.4 Pa in "1989" seasons for Triton 891 !! converges toward 1.1 Pa in "2015" seasons for Pluto 892 if (convergeps) then 893 if (triton) then 894 ! 1989 declination 895 if (declin*180./pi.gt.-46..and.declin*180./pi.lt.-45. & 896 .and.zday.gt.saveday+1000. & 897 .and.declin.lt.savedeclin) then 898 call globalaverage2d(ngrid,pplev(:,1),globave) 899 if (globave.gt.1.5) then 900 adjust=adjust+0.005 901 else if (globave.lt.1.3) then 902 adjust=adjust-0.005 903 endif 904 saveday=zday 905 endif 906 else 907 ! Pluto : 2015 declination current epoch 908 if (declin*180./pi.gt.50.and.declin*180./pi.lt.51. & 909 .and.zday.gt.saveday+10000. & 910 .and.declin.gt.savedeclin) then 911 call globalaverage2d(ngrid,pplev(:,1),globave) 912 if (globave.gt.1.2) then 913 adjust=adjust+0.005 914 else if (globave.lt.1.) then 915 adjust=adjust-0.005 916 endif 917 saveday=zday 918 endif 919 endif 920 end if 921 end if ! if ngrid ne 1 922 923 call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat, & 924 capcal,adjust,dist_star,albedo,emis,flusurfold,ptimestep, & 925 zls) 926 ! do i=2,ngrid 927 ! albedo(i,:) = albedo(1,:) 928 ! enddo 878 929 879 930 if (callrad) then … … 920 971 921 972 ! standard callcorrk 922 ! clearsky=.false.923 ! if (oldplutocorrk) then924 ! call callcorrk_pluto(icount,ngrid,nlayer,pq,nq,qsurf,&925 ! albedo,emis,mu0,pplev,pplay,pt,&926 ! tsurf,fract,dist_star,aerosol,&927 ! zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,&928 ! fluxabs_sw,fluxtop_dn,reffrad,tau_col,ptime,pday,&929 ! firstcall,lastcall,zzlay)930 !else931 call callcorrk(ngrid,nlayer,pq,nq,qsurf, 973 if (oldplutocorrk) then 974 call callcorrk_pluto(icount,ngrid,nlayer,pq,nq,qsurf, & 975 albedo(:,1),emis,mu0,pplev,pplay,pt, & 976 zzlay,tsurf,fract,dist_star,aerosol, & 977 zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, & 978 fluxabs_sw,fluxtop_dn,reffrad,tau_col,ptime,pday, & 979 cloudfrac,totcloudfrac,.false., & 980 firstcall,lastcall) 981 else 982 call callcorrk(ngrid,nlayer,pq,nq,qsurf, & 932 983 albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & 933 984 zzlay,tsurf,fract,dist_star,aerosol,muvar, & … … 938 989 tau_col,cloudfrac,totcloudfrac, & 939 990 .false.,firstcall,lastcall) 940 !endif ! oldplutocorrk991 endif ! oldplutocorrk 941 992 !GG (feb2021): Option to "artificially" decrease the raditive time scale in 942 993 !the deep atmosphere press > 0.1 bar. Suggested by MT. … … 956 1007 957 1008 ! AF24: removed CLFvarying Option 958 959 ! if(ok_slab_ocean) then960 ! tsurf(:)=tsurf2(:)961 ! endif962 963 1009 964 1010 ! Radiative flux from the sky absorbed by the surface (W.m-2). … … 1052 1098 1053 1099 if (oldplutovdifc) then 1054 zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)1055 1100 zdum1(:,:) = 0 1056 1101 zdum2(:,:) = 0 … … 1066 1111 zdqdif,zdqsdif,qsat_ch4,qsat_ch4_l1) !,zq1temp_ch4,qsat_ch4) 1067 1112 1068 pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)1069 pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)1070 pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)1071 1072 1113 zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only 1073 zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)1074 1114 1075 1115 bcond=1./tcond1p4Pa 1076 1116 acond=r/lw_n2 1077 1078 if (tracer) then1079 pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)1080 dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)1081 end if ! of if (tracer)1082 1117 1083 1118 !------------------------------------------------------------------ … … 1379 1414 zdqphot_prec(:,:)=0. 1380 1415 zdqphot_ch4(:,:)=0. 1416 zdqhaze(:,:,:)=0 1381 1417 ! Forcing to a fixed haze profile if haze_proffix 1382 1418 if (haze_proffix.and.i_haze.gt.0.) then 1383 call haze_prof(ngrid,nlayer,zzlay,pplay,pt, &1384 1385 zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze)) &1386 1419 call haze_prof(ngrid,nlayer,zzlay,pplay,pt, & 1420 reffrad,profmmr) 1421 zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze)) & 1422 /ptimestep 1387 1423 else 1388 call hazecloud(ngrid,nlayer,nq,ptimestep, &1389 pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,zdqhaze, &1390 zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)1424 call hazecloud(ngrid,nlayer,nq,ptimestep, & 1425 pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,zdqhaze, & 1426 zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin) 1391 1427 endif 1392 1428 … … 1589 1625 ! todo: should be placed in haze and use tendency of n2 instead of flusurf 1590 1626 IF (source_haze) THEN 1627 write(*,*) "Source haze not supported yet." 1628 stop 1591 1629 ! call hazesource(ngrid,nlayer,nq,ptimestep, & 1592 1630 ! pplev,flusurf,mu0,zdq_source) -
trunk/LMDZ.PLUTO/libf/phypluto/radii_mod.F90
r3195 r3329 28 28 use ioipsl_getin_p_mod, only: getin_p 29 29 use radinc_h, only: naerkind 30 use datafile_mod, only: hazerad_file 30 31 use aerosol_mod, only: iaero_haze, i_haze, & 31 32 iaero_generic, i_rgcs_ice … … 50 51 !reffrad(1:ngrid,1:nlayer,iaer) = 2.e-6 ! haze 51 52 nueffrad(1:ngrid,1:nlayer,iaer) = 0.02 ! haze 52 print*, 'TB22 Hello2',radius(i_haze)*nmono**(1./3.)53 53 endif 54 54 enddo … … 100 100 IF (firstcall) then 101 101 firstcall=.false. 102 file_path=trim(datadir)//'/haze_prop/ hazerad.txt'102 file_path=trim(datadir)//'/haze_prop/'//hazerad_file 103 103 open(223,file=file_path,form='formatted') 104 104 do ifine=1,Nfine -
trunk/LMDZ.PLUTO/libf/phypluto/suaer_corrk.F90
r3195 r3329 10 10 use radinc_h, only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind 11 11 use radcommon_h, only: blamv,blami,lamrefir,lamrefvis 12 use datafile_mod, only: datadir, aerdir, hazeprop 12 use datafile_mod, only: datadir, aerdir, hazeprop_file 13 13 14 14 ! outputs … … 142 142 ! Only one table of optical properties : 143 143 write(*,*)'Suaer haze optical properties, using: ', & 144 TRIM(hazeprop )144 TRIM(hazeprop_file) 145 145 146 146 ! visible 147 file_id(iaer,1)=TRIM(hazeprop )147 file_id(iaer,1)=TRIM(hazeprop_file) 148 148 ! infrared 149 149 file_id(iaer,2)=file_id(iaer,1) … … 155 155 enddo 156 156 157 IF (naerkind .gt. 1) THEN158 print*,'naerkind = ',naerkind159 print*,'but we only have data for 1 type, exiting.'160 call abort161 ENDIF157 ! IF (naerkind .gt. 1) THEN ! AF24: ? 158 ! print*,'naerkind = ',naerkind 159 ! print*,'but we only have data for 1 type, exiting.' 160 ! call abort 161 ! ENDIF 162 162 163 163 !------------------------------------------------------------------ … … 314 314 315 315 jfile = jfile+1 316 IF ((idomain.NE.iaero_haze).OR.(iaer.NE.iaero_haze)) THEN 316 print *, idomain, " ", iaer, " ", iaero_haze 317 IF ((idomain.NE.1).OR.(iaer.NE.1).OR.(jfile.EQ.6)) THEN 318 ! AF24: what does it mean? :-O 317 319 endwhile = .true. 318 320 ENDIF -
trunk/LMDZ.PLUTO/libf/phypluto/surfprop.F90
r3275 r3329 4 4 5 5 ! use comgeomfi_h, only: 6 use radinc_h, only: L_NSPECTV, L_NSPECTI 6 7 use comcstfi_mod, only: pi 7 8 use comsoil_h, only: nsoilmx,mlayer,volcapa,inertiedat … … 64 65 REAL,INTENT(IN) :: adjust 65 66 REAL,INTENT(IN) :: dist 66 REAL,INTENT(OUT) :: albedo(ngrid )67 REAL,INTENT(OUT) :: albedo(ngrid,L_NSPECTV) 67 68 REAL,INTENT(OUT) :: emis(ngrid) 68 69 REAL,INTENT(OUT) :: tidat(ngrid,nsoilmx) … … 169 170 170 171 CASE (0) ! fixed albedo 171 albedo(ig )=min(max(alb_n2b+adjust,0.),1.)172 albedo(ig,:)=min(max(alb_n2b+adjust,0.),1.) 172 173 173 174 CASE (1) ! Albedo decreases with thickness 174 albedo(ig )=(alb_n2b-deltab)/(1.-10000.)*qsurf(ig,igcm_n2)+alb_n2b175 albedo(ig )=min(max(alb_n2b-deltab,albedo(ig)),alb_n2b) ! should not be higher than alb_n2b and not lower than alb_n2b-deltab175 albedo(ig,:)=(alb_n2b-deltab)/(1.-10000.)*qsurf(ig,igcm_n2)+alb_n2b 176 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 176 177 CASE (2) ! Special Sputnik differences of albedo 177 178 if (qsurf(ig,igcm_n2).ge.1.e4.and.(longitude(ig)*180./pi.le.-161.or.longitude(ig)*180./pi.ge.146)) then … … 179 180 ((longitude(ig)*180./pi.gt.140.).and. & 180 181 (longitude(ig)*180./pi.lt.165.)) ) then 181 albedo(ig )=alb_n2b-deltab182 albedo(ig,:)=alb_n2b-deltab 182 183 else 183 albedo(ig )=alb_n2b+deltab184 albedo(ig,:)=alb_n2b+deltab 184 185 endif 185 186 else 186 albedo(ig )=alb_n2b187 albedo(ig,:)=alb_n2b 187 188 endif 188 189 189 190 CASE (3) ! Albedo increases (delta neg) or decreases (delta pos) with sublimation rates 190 albedo(ig )=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig)191 albedo(ig )=min(max(alb_n2b-deltab,albedo(ig)),alb_n2b+deltab)191 albedo(ig,:)=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig,:) 192 albedo(ig,:)=min(max(alb_n2b-deltab,albedo(ig,:)),alb_n2b+deltab) 192 193 !! Triton: 193 !albedo(ig )=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig)194 !albedo(ig )=min(max(alb_n2b-deltab,albedo(ig)),alb_n2b+deltab)194 !albedo(ig,:)=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig,:) 195 !albedo(ig,:)=min(max(alb_n2b-deltab,albedo(ig,:)),alb_n2b+deltab) 195 196 196 197 CASE (4) ! Albedo Difference in N/S 197 198 if (latitude(ig)*180./pi.ge.0.) then 198 albedo(ig )=min(max(alb_n2b-deltab+adjust,0.),1.)199 albedo(ig,:)=min(max(alb_n2b-deltab+adjust,0.),1.) 199 200 else 200 albedo(ig )=min(max(alb_n2b+deltab+adjust,0.),1.)201 albedo(ig,:)=min(max(alb_n2b+deltab+adjust,0.),1.) 201 202 endif 202 203 … … 210 211 ! .or. ((latitude(ig)*180./pi.le.6.).and.(latitude(ig)*180./pi.ge.2.)) & 211 212 ! ) then 212 ! albedo(ig )=alb_n2b-deltab213 ! albedo(ig,:)=alb_n2b-deltab 213 214 214 215 !! Patches GCM … … 218 219 .or. ((latitude(ig)*180./pi.lt.30.5).and.(latitude(ig)*180./pi.gt.29.5).and. & 219 220 (longitude(ig)*180./pi.lt.169.).and.(longitude(ig)*180./pi.gt.168.))) then 220 albedo(ig )=alb_n2b-deltab221 albedo(ig,:)=alb_n2b-deltab 221 222 else if (((latitude(ig)*180./pi.lt.30.5).and.(latitude(ig)*180./pi.gt.29.).and. & 222 223 (longitude(ig)*180./pi.lt.165.5).and.(longitude(ig)*180./pi.gt.164.5)) & 223 224 .or. ((latitude(ig)*180./pi.lt.33.).and.(latitude(ig)*180./pi.gt.32.).and. & 224 225 (longitude(ig)*180./pi.lt.169.).and.(longitude(ig)*180./pi.gt.168.))) then 225 albedo(ig )=alb_n2b+deltab226 albedo(ig,:)=alb_n2b+deltab 226 227 else 227 albedo(ig )=alb_n2b228 albedo(ig,:)=alb_n2b 228 229 endif 229 230 else 230 albedo(ig )=alb_n2b231 albedo(ig,:)=alb_n2b 231 232 endif 232 233 233 234 CASE (7) ! Albedo from albedodat and adjusted emissivity 234 albedo(ig )=albedodat(ig)235 albedo(ig,:)=albedodat(ig) 235 236 ! adjust emissivity if convergeps = true 236 237 emis(ig)=min(max(emis(ig)+adjust,0.),1.) … … 247 248 248 249 else if (ice_case.eq.2) then 249 albedo(ig )=alb_co250 albedo(ig,:)=alb_co 250 251 emis(ig)=emis_co 251 252 … … 261 262 emis(ig)=emis_ch4 262 263 if (latitude(ig)*180./pi.le.metlateq.and.latitude(ig)*180./pi.gt.-metlateq) then 263 albedo(ig )=alb_ch4_eq264 albedo(ig,:)=alb_ch4_eq 264 265 else 265 albedo(ig )=alb_ch4266 albedo(ig,:)=alb_ch4 266 267 endif 267 268 … … 269 270 emis(ig)=emis_ch4 270 271 if (latitude(ig)*180./pi.le.metlateq.and.latitude(ig)*180./pi.gt.-metlateq) then 271 albedo(ig )=alb_ch4_eq272 albedo(ig,:)=alb_ch4_eq 272 273 else if (latitude(ig)*180./pi.le.-metlateq) then 273 albedo(ig )=alb_ch4_s274 albedo(ig,:)=alb_ch4_s 274 275 else 275 albedo(ig )=alb_ch4276 albedo(ig,:)=alb_ch4 276 277 endif 277 278 … … 279 280 280 281 emis(ig)=emis_ch4 281 albedo(ig )=alb_ch4282 albedo(ig,:)=alb_ch4 282 283 283 284 if (latitude(ig)*180./pi.le.metlateq.and.latitude(ig)*180./pi.gt.-metlateq) then 284 albedo(ig )=alb_ch4_eq285 albedo(ig,:)=alb_ch4_eq 285 286 endif 286 287 … … 291 292 SELECT CASE (feedback_met) 292 293 CASE (0) ! Default (linear from alb_ch4_eq) 293 albedo(ig )=(1.-alb_ch4_eq)/0.002*qsurf(ig,igcm_ch4_ice)+alb_ch4_eq294 albedo(ig,:)=(1.-alb_ch4_eq)/0.002*qsurf(ig,igcm_ch4_ice)+alb_ch4_eq 294 295 !emis(ig)=(1.-emis_ch4)/0.002*qsurf(ig,igcm_ch4_ice)+emis_ch4 295 if (albedo(ig ).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb296 if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb 296 297 !if (emis(ig).gt.1.) emis(ig)=1. 297 298 CASE (1) ! Hyperbolic tangent old 298 albedo(ig )=0.6*0.5*(1.+tanh(6.*(0.001+qsurf(ig,igcm_ch4_ice))/0.5))+0.3 !299 albedo(ig,:)=0.6*0.5*(1.+tanh(6.*(0.001+qsurf(ig,igcm_ch4_ice))/0.5))+0.3 ! 299 300 !emis(ig)=0.2*0.5*(1.+tanh(6.*(0.001+qsurf(ig,igcm_ch4_ice))/0.5))+0.8 ! 300 if (albedo(ig ).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb301 if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb 301 302 !if (emis(ig).gt.1.) emis(ig)=1. 302 303 CASE (2) ! hyperbolic tangent old 303 albedo(ig )=0.5*0.6*(1.+tanh(1000.*(qsurf(ig,igcm_ch4_ice)-0.002)))+0.3 !304 albedo(ig,:)=0.5*0.6*(1.+tanh(1000.*(qsurf(ig,igcm_ch4_ice)-0.002)))+0.3 ! 304 305 !emis(ig)=0.2*0.5*(1.+tanh(1000.*(qsurf(ig,igcm_ch4_ice)-0.002)))+0.8 ! 305 if (albedo(ig ).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb306 if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb 306 307 !if (emis(ig).gt.1.) emis(ig)=1. 307 308 CASE (3) ! hyperbolic tangent equation with parameters 308 albedo(ig )=aa*(-1+tanh(fdch4_ampl*qsurf(ig,igcm_ch4_ice)))+bb309 if (albedo(ig ).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb309 albedo(ig,:)=aa*(-1+tanh(fdch4_ampl*qsurf(ig,igcm_ch4_ice)))+bb 310 if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb 310 311 CASE DEFAULT 311 312 write(*,*) 'STOP surfprop.F90: fd wrong' … … 319 320 emis(ig)=emis_ch4 320 321 if (latitude(ig)*180./pi.le.metlateq.and.latitude(ig)*180./pi.gt.-metlateq) then 321 albedo(ig )=alb_ch4_eq ! Pure methane ice322 albedo(ig,:)=alb_ch4_eq ! Pure methane ice 322 323 else if (latitude(ig)*180./pi.le.-metlateq) then 323 albedo(ig )=alb_ch4_s ! Pure methane ice324 albedo(ig,:)=alb_ch4_s ! Pure methane ice 324 325 if (pls.le.metls2.and.pls.gt.metls1) then 325 albedo(ig )=alb_ch4_s+deltab ! Also used for N2, careful326 albedo(ig,:)=alb_ch4_s+deltab ! Also used for N2, careful 326 327 endif 327 328 else 328 albedo(ig )=alb_ch4 ! Pure methane ice329 albedo(ig,:)=alb_ch4 ! Pure methane ice 329 330 endif 330 331 331 332 CASE (4) ! Albedo from albedodat 332 333 emis(ig)=emis_ch4 333 albedo(ig )=albedodat(ig)334 albedo(ig,:)=albedodat(ig) 334 335 335 336 CASE DEFAULT … … 349 350 350 351 if (latitude(ig)*180./pi.le.tholateq.and.latitude(ig)*180./pi.gt.-tholateq) then 351 albedo(ig )=alb_tho_eq352 albedo(ig,:)=alb_tho_eq 352 353 emis(ig)=emis_tho_eq 353 354 else 354 albedo(ig )=alb_tho ! default = 0.1355 albedo(ig,:)=alb_tho ! default = 0.1 355 356 emis(ig)=emis_tho 356 357 endif … … 359 360 360 361 if (latitude(ig)*180./pi.le.tholateq.and.latitude(ig)*180./pi.gt.-tholateq) then 361 albedo(ig )=alb_tho_eq362 albedo(ig,:)=alb_tho_eq 362 363 emis(ig)=emis_tho_eq 363 364 else 364 albedo(ig )=alb_tho ! default = 0.1365 albedo(ig,:)=alb_tho ! default = 0.1 365 366 emis(ig)=emis_tho 366 367 endif … … 368 369 if (latitude(ig)*180./pi.le.tholatn.and.latitude(ig)*180./pi.ge.tholats & 369 370 .and.longitude(ig)*180./pi.ge.tholone.and.longitude(ig)*180./pi.le.tholonw) then 370 albedo(ig )=alb_tho_spe371 albedo(ig,:)=alb_tho_spe 371 372 emis(ig)=emis_tho_spe 372 373 endif … … 374 375 CASE (2) ! Depends on the albedo map 375 376 376 albedo(ig )=albedodat(ig)377 if (albedo(ig ).gt.alb_ch4) then377 albedo(ig,:)=albedodat(ig) 378 if (albedo(ig,1).gt.alb_ch4) then 378 379 emis(ig)=emis_ch4 379 380 else … … 391 392 392 393 if (specalb) then 393 albedo(ig )=albedodat(ig) ! Specific ground properties394 albedo(ig,:)=albedodat(ig) ! Specific ground properties 394 395 !if (albedodat(ig).lt.0.25) then 395 ! albedo(ig )=alb_tho_eq396 ! albedo(ig,:)=alb_tho_eq 396 397 !else if (albedodat(ig).lt.0.40) then 397 ! albedo(ig )=0.35398 ! albedo(ig,:)=0.35 398 399 !else if (albedodat(ig).lt.0.70) then 399 ! albedo(ig )=0.51400 ! albedo(ig,:)=0.51 400 401 !endif 401 402 endif
Note: See TracChangeset
for help on using the changeset viewer.