Changeset 1974 for trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F
- Timestamp:
- Jul 18, 2018, 4:48:34 PM (6 years ago)
- File:
-
- 1 edited
-
trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F
r1861 r1974 7 7 SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls, 8 8 & pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad, 9 & nueffrad,QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d, 9 & QREFvis3d,QREFir3d,omegaREFir3d, 10 & totstormfract,clearatm,dsords, 10 11 & clearsky,totcloudfrac) 11 12 12 13 ! to use 'getin' 13 14 USE ioipsl_getincom, only: getin 14 15 use tracer_mod, only: noms, igcm_h2o_ice, igcm_dust_mass, 15 16 & igcm_dust_submicron, rho_dust, rho_ice, 16 & nqdust 17 & nqdust, igcm_stormdust_mass 17 18 use geometry_mod, only: latitude ! grid point latitudes (rad) 18 19 use comgeomfi_h, only: sinlat ! sines of grid point latitudes … … 26 27 & iaerdust,tauvis, 27 28 & iaer_dust_conrath,iaer_dust_doubleq, 28 & iaer_dust_submicron,iaer_h2o_ice 29 & iaer_dust_submicron,iaer_h2o_ice, 30 & iaer_stormdust_doubleq 31 USE calcstormfract_mod 29 32 IMPLICIT NONE 30 33 c======================================================================= … … 56 59 c QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients 57 60 c QREFir3d(ngrid,nlayer,naerkind) / at reference wavelengths; 58 c omegaREFvis3d(ngrid,nlayer,naerkind) \ 3d single scat. albedo59 61 c omegaREFir3d(ngrid,nlayer,naerkind) / at reference wavelengths; 60 62 c … … 65 67 c aerosol aerosol(ig,l,1) is the dust optical 66 68 c depth in layer l, grid point ig 67 69 c taualldust CW17 total opacity for all dust scatterer stormdust included 68 70 c 69 71 c======================================================================= 70 #include "callkeys.h"72 include "callkeys.h" 71 73 72 74 c----------------------------------------------------------------------- … … 77 79 c Input/Output 78 80 c ------------ 79 INTEGER ngrid,nlayer,nq80 81 REAL ls,zday,expfactor82 REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)83 REAL pq(ngrid,nlayer,nq)84 REAL tauref(ngrid),tau(ngrid,naerkind)85 REAL aerosol(ngrid,nlayer,naerkind)86 REAL dsodust(ngrid,nlayer)87 REAL reffrad(ngrid,nlayer,naerkind)88 REAL nueffrad(ngrid,nlayer,naerkind)89 REAL QREFvis3d(ngrid,nlayer,naerkind)90 REAL QREFir3d(ngrid,nlayer,naerkind)91 REAL omegaREFvis3d(ngrid,nlayer,naerkind)92 REAL omegaREFir3d(ngrid,nlayer,naerkind)93 real,intent(in) :: totcloudfrac(ngrid) ! total cloud fraction94 logical,intent(in) :: clearsky ! true for part without clouds,95 ! false for part with clouds (total or sub-grid clouds)96 real :: CLFtot ! total cloud fraction81 INTEGER, INTENT(IN) :: ngrid,nlayer,nq 82 REAL, INTENT(IN) :: ls,zday 83 REAL, INTENT(IN) :: pplev(ngrid,nlayer+1),pplay(ngrid,nlayer) 84 REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) 85 REAL, INTENT(OUT) :: tauref(ngrid) 86 REAL, INTENT(OUT) :: tau(ngrid,naerkind) 87 REAL, INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) 88 REAL, INTENT(OUT) :: dsodust(ngrid,nlayer) 89 REAL, INTENT(OUT) :: dsords(ngrid,nlayer) !dso of stormdust 90 REAL, INTENT(INOUT) :: reffrad(ngrid,nlayer,naerkind) 91 REAL, INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) 92 REAL, INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind) 93 REAL, INTENT(IN) :: omegaREFir3d(ngrid,nlayer,naerkind) 94 LOGICAL, INTENT(IN) :: clearatm 95 REAL, INTENT(IN) :: totstormfract(ngrid) 96 REAL, INTENT(OUT) :: tauscaling(ngrid) ! Scaling factor for qdust and Ndust 97 REAL,INTENT(IN) :: totcloudfrac(ngrid) ! total cloud fraction 98 LOGICAL,INTENT(IN) :: clearsky ! true for part without clouds,false for part with clouds (total or sub-grid clouds) 97 99 c 98 100 c Local variables : 99 101 c ----------------- 102 REAL CLFtot ! total cloud fraction 103 real expfactor 100 104 INTEGER l,ig,iq,i,j 101 105 INTEGER iaer ! Aerosol index … … 110 114 c for dust or water ice particles in the radiative transfer 111 115 c (see callradite.F for more information). 112 REAL taudusttmp(ngrid)! Temporary dust opacity 113 ! used before scaling 114 REAL tauscaling(ngrid) ! Scaling factor for qdust and Ndust 116 REAL taudusttmp(ngrid)! Temporary dust opacity used before scaling 117 REAL taubackdusttmp(ngrid)! Temporary background dust opacity used before scaling 118 REAL taualldust(ngrid)! dust opacity all dust 119 REAL taudust(ngrid)! dust opacity dust doubleq 120 REAL taustormdust(ngrid)! dust opacity stormdust doubleq 121 REAL taustormdusttmp(ngrid)! dust opacity stormdust doubleq before tauscaling 115 122 REAL taudustvis(ngrid) ! Dust opacity after scaling 116 123 REAL taudusttes(ngrid) ! Dust opacity at IR ref. wav. as … … 157 164 158 165 tau(1:ngrid,1:naerkind)=0 166 dsords(:,:)=0. !CW17: initialize dsords 167 dsodust(:,:)=0. 159 168 160 169 ! identify tracers … … 166 175 DO iaer=1,naerkind 167 176 txt=name_iaer(iaer) 168 IF (txt(1:4).eq."dust") THEN 177 ! CW17: choice tauscaling for stormdust or not 178 IF ((txt(1:4).eq."dust").OR.(txt(1:5).eq."storm")) THEN 169 179 naerdust=naerdust+1 170 180 iaerdust(naerdust)=iaer … … 213 223 call getin("tauvis",tauvis) 214 224 215 IF (freedust ) THEN225 IF (freedust.or.rdstorm) THEN ! if rdstorm no need to held opacity constant at the first levels 216 226 cstdustlevel = 1 217 227 ELSE 218 cstdustlevel = cstdustlevel0 228 cstdustlevel = cstdustlevel0 !Opacity in the first levels is held constant to 229 !avoid unrealistic values due to constant lifting 219 230 ENDIF 220 231 … … 374 385 & ( pplev(ig,l) - pplev(ig,l+1) ) / g 375 386 ! DENSITY SCALED OPACITY IN INFRARED: 376 dsodust(ig,l) =377 & ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /378 & ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *379 & pq(ig,cstdustlevel,igcm_dust_mass)387 ! dsodust(ig,l) = 388 ! & ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) / 389 ! & ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) * 390 ! & pq(ig,cstdustlevel,igcm_dust_mass) 380 391 ENDDO 381 392 ELSE … … 387 398 & ( pplev(ig,l) - pplev(ig,l+1) ) / g 388 399 ! DENSITY SCALED OPACITY IN INFRARED: 389 dsodust(ig,l) =390 & ( 0.75 * QREFir3d(ig,l,iaer) /391 & ( rho_dust * reffrad(ig,l,iaer) ) ) *392 & pq(ig,l,igcm_dust_mass)400 ! dsodust(ig,l) = 401 ! & ( 0.75 * QREFir3d(ig,l,iaer) / 402 ! & ( rho_dust * reffrad(ig,l,iaer) ) ) * 403 ! & pq(ig,l,igcm_dust_mass) 393 404 ENDDO 394 405 ENDIF … … 473 484 ENDIF ! end (clearsky) 474 485 475 c 3. Outputs -- Now done in physiq.F 476 ! IF (ngrid.NE.1) THEN 477 ! CALL WRITEDIAGFI(ngrid,'tauVIS','tauext VIS refwvl', 478 ! & ' ',2,taucloudvis) 479 ! CALL WRITEDIAGFI(ngrid,'tauTES','tauabs IR refwvl', 480 ! & ' ',2,taucloudtes) 481 ! IF (callstats) THEN 482 ! CALL wstats(ngrid,'tauVIS','tauext VIS refwvl', 483 ! & ' ',2,taucloudvis) 484 ! CALL wstats(ngrid,'tauTES','tauabs IR refwvl', 485 ! & ' ',2,taucloudtes) 486 ! ENDIF 487 ! ELSE 488 ! CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU') 489 ! ENDIF 486 c================================================================== 487 CASE("stormdust_doubleq") aerkind ! CW17 : Two-moment scheme for 488 c stormdust (transport of mass and number mixing ratio) 489 c================================================================== 490 c aerosol is calculated twice : once within the dust storm (clearatm=false) 491 c and once in the part of the mesh without dust storm (clearatm=true) 492 aerosol(1:ngrid,1:nlayer,iaer) = 0. 493 IF (clearatm) THEN ! considering part of the mesh without storm 494 aerosol(1:ngrid,1:nlayer,iaer)=1.e-25 495 ELSE ! part of the mesh with concentred dust storm 496 DO l=1,nlayer 497 IF (l.LE.cstdustlevel) THEN 498 c Opacity in the first levels is held constant to 499 c avoid unrealistic values due to constant lifting: 500 DO ig=1,ngrid 501 aerosol(ig,l,iaer) = 502 & ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) / 503 & ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) * 504 & pq(ig,cstdustlevel,igcm_stormdust_mass) * 505 & ( pplev(ig,l) - pplev(ig,l+1) ) / g 506 ENDDO 507 ELSE 508 DO ig=1,ngrid 509 aerosol(ig,l,iaer) = 510 & ( 0.75 * QREFvis3d(ig,l,iaer) / 511 & ( rho_dust * reffrad(ig,l,iaer) ) ) * 512 & pq(ig,l,igcm_stormdust_mass) * 513 & ( pplev(ig,l) - pplev(ig,l+1) ) / g 514 ENDDO 515 ENDIF 516 ENDDO 517 ENDIF 490 518 c================================================================== 491 519 END SELECT aerkind … … 498 526 c is originally scaled to an equivalent odpref Pa pressure surface. 499 527 c ----------------------------------------------------------------- 528 500 529 501 530 #ifdef DUSTSTORM … … 640 669 IF (freedust) THEN 641 670 tauscaling(:) = 1. 642 671 c opacity obtained with stormdust 672 IF (rdstorm) THEN 673 taustormdusttmp(1:ngrid)=0. 674 DO l=1,nlayer 675 DO ig=1,ngrid 676 taustormdusttmp(ig) = taustormdusttmp(ig)+ 677 & aerosol(ig,l,iaerdust(2)) 678 ENDDO 679 ENDDO 680 !opacity obtained with background dust only 681 taubackdusttmp(1:ngrid)=0. 682 DO l=1,nlayer 683 DO ig=1,ngrid 684 taubackdusttmp(ig) = taubackdusttmp(ig)+ 685 & aerosol(ig,l,iaerdust(1)) 686 ENDDO 687 ENDDO 688 ENDIF !rdsstorm 643 689 ELSE 644 690 c Temporary scaling factor … … 669 715 & aerosol(ig,l,iaerdust(iaer))* tauscaling(ig)) 670 716 ENDDO 671 ENDDO672 ENDDO673 674 DO l=1,nlayer675 DO ig=1,ngrid676 dsodust(ig,l) = max(1E-20,dsodust(ig,l)* tauscaling(ig))677 717 ENDDO 678 718 ENDDO … … 700 740 tauref(:) = tauref(:) * odpref / pplev(:,1) 701 741 ENDIF 702 703 c output for debug 704 c IF (ngrid.NE.1) THEN 705 c CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust', 706 c & '#',2,taudusttmp) 707 c CALL WRITEDIAGFI(ngrid,'tausca','tauscaling', 708 c & '#',2,tauscaling) 709 c ELSE 710 c CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust', 711 c & '#',0,taudusttmp) 712 c CALL WRITEDIAGFI(ngrid,'tausca','tauscaling', 713 c & '#',0,tauscaling) 714 c ENDIF 742 715 743 c ----------------------------------------------------------------- 716 744 c Column integrated visible optical depth in each point … … 724 752 ENDDO 725 753 754 c for diagnostics: opacity for all dust scatterers stormdust included 755 taualldust(1:ngrid)=0. 756 DO iaer=1,naerdust 757 DO l=1,nlayer 758 DO ig=1,ngrid 759 taualldust(ig) = taualldust(ig) + 760 & aerosol(ig,l,iaerdust(iaer)) 761 ENDDO 762 ENDDO 763 ENDDO 764 765 IF (rdstorm) THEN 766 767 c for diagnostics: opacity for dust in background only 768 taudust(1:ngrid)=0. 769 DO l=1,nlayer 770 DO ig=1,ngrid 771 taudust(ig) = taudust(ig) + 772 & aerosol(ig,l,iaer_dust_doubleq) 773 ENDDO 774 ENDDO 775 776 c for diagnostics: opacity for dust in storm only 777 taustormdust(1:ngrid)=0. 778 DO l=1,nlayer 779 DO ig=1,ngrid 780 taustormdust(ig) = taustormdust(ig) + 781 & aerosol(ig,l,iaer_stormdust_doubleq) 782 ENDDO 783 ENDDO 784 785 ENDIF 786 787 726 788 #ifdef DUSTSTORM 727 789 IF (firstcall) THEN … … 733 795 c Density scaled opacity and column opacity output 734 796 c ----------------------------------------------------------------- 735 c dsodust(1:ngrid,1:nlayer) = 0. 736 c DO iaer=1,naerdust 737 c DO l=1,nlayer 738 c DO ig=1,ngrid 739 c dsodust(ig,l) = dsodust(ig,l) + 740 c & aerosol(ig,l,iaerdust(iaer)) * g / 741 c & (pplev(ig,l) - pplev(ig,l+1)) 742 c ENDDO 743 c ENDDO 744 c IF (ngrid.NE.1) THEN 745 c write(txt2,'(i1.1)') iaer 746 c call WRITEDIAGFI(ngrid,'taudust'//txt2, 747 c & 'Dust col opacity', 748 c & ' ',2,tau(1,iaerdust(iaer))) 749 c IF (callstats) THEN 750 c CALL wstats(ngrid,'taudust'//txt2, 751 c & 'Dust col opacity', 752 c & ' ',2,tau(1,iaerdust(iaer))) 753 c ENDIF 754 c ENDIF 755 c ENDDO 756 757 c IF (ngrid.NE.1) THEN 758 c CALL WRITEDIAGFI(ngrid,'dsodust','tau*g/dp', 759 c & 'm2.kg-1',3,dsodust) 760 c IF (callstats) THEN 761 c CALL wstats(ngrid,'dsodust', 762 c & 'tau*g/dp', 763 c & 'm2.kg-1',3,dsodust) 764 c ENDIF 765 c ELSE 766 c CALL WRITEDIAGFI(ngrid,"dsodust","dsodust","m2.kg-1",1, 767 c & dsodust) 768 c ENDIF ! of IF (ngrid.NE.1) 769 c ----------------------------------------------------------------- 797 IF (rdstorm) then 798 DO l=1,nlayer 799 IF (l.LE.cstdustlevel) THEN 800 DO ig=1,ngrid 801 dsodust(ig,l)=dsodust(ig,l) + 802 & aerosol(ig,l,iaer_dust_doubleq) * g / 803 & (pplev(ig,l) - pplev(ig,l+1)) 804 805 dsords(ig,l) = dsords(ig,l) + 806 & aerosol(ig,l,iaer_stormdust_doubleq)* g/ 807 & (pplev(ig,l) - pplev(ig,l+1)) 808 ENDDO 809 ELSE 810 DO ig=1,ngrid 811 dsodust(ig,l) =dsodust(ig,l) + 812 & aerosol(ig,l,iaer_dust_doubleq) * g / 813 & (pplev(ig,l) - pplev(ig,l+1)) 814 dsords(ig,l) = dsords(ig,l) + 815 & aerosol(ig,l,iaer_stormdust_doubleq)* g/ 816 & (pplev(ig,l) - pplev(ig,l+1)) 817 ENDDO 818 ENDIF 819 ENDDO 820 ENDIF 821 822 c ----------------------------------------------------------------- 823 c ----------------------------------------------------------------- 824 c aerosol/X for stormdust to prepare calculation of radiative transfer 825 c ----------------------------------------------------------------- 826 if (rdstorm) then 827 DO l=1,nlayer 828 DO ig=1,ngrid 829 aerosol(ig,l,iaer_stormdust_doubleq) = 830 & aerosol(ig,l,iaer_stormdust_doubleq)/totstormfract(ig) 831 ENDDO 832 ENDDO 833 endif 834 770 835 771 836 END SUBROUTINE aeropacity
Note: See TracChangeset
for help on using the changeset viewer.
