!$Id: photo_interp.F90 104 2008-12-23 10:28:51Z acosce $ !! ========================================================================= !! INCA - INteraction with Chemistry and Aerosols !! !! Copyright Laboratoire des Sciences du Climat et de l'Environnement (LSCE) !! Unite mixte CEA-CNRS-UVSQ !! !! Contributors to this INCA subroutine: !! !! Stacy Walters, NCAR, stacy@ucar.edu !! !! Anne Cozic, LSCE, anne.cozic@cea.fr !! Yann Meurdesoif, LSCE, yann.meurdesoif@cea.fr !! !! This software is a computer program whose purpose is to simulate the !! atmospheric gas phase and aerosol composition. The model is designed to be !! used within a transport model or a general circulation model. This version !! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts !! for emissions, transport (resolved and sub-grid scale), photochemical !! transformations, and scavenging (dry deposition and washout) of chemical !! species and aerosols interactively in the GCM. Several versions of the INCA !! model are currently used depending on the envisaged applications with the !! chemistry-climate model. !! !! This software is governed by the CeCILL license under French law and !! abiding by the rules of distribution of free software. You can use, !! modify and/ or redistribute the software under the terms of the CeCILL !! license as circulated by CEA, CNRS and INRIA at the following URL !! "http://www.cecill.info". !! !! As a counterpart to the access to the source code and rights to copy, !! modify and redistribute granted by the license, users are provided only !! with a limited warranty and the software's author, the holder of the !! economic rights, and the successive licensors have only limited !! liability. !! !! In this respect, the user's attention is drawn to the risks associated !! with loading, using, modifying and/or developing or reproducing the !! software by the user in light of its specific status of free software, !! that may mean that it is complicated to manipulate, and that also !! therefore means that it is reserved for developers and experienced !! professionals having in-depth computer knowledge. Users are therefore !! encouraged to load and test the software's suitability as regards their !! requirements in conditions enabling the security of their systems and/or !! data to be ensured and, more generally, to use and operate it in the !! same conditions as regards security. !! !! The fact that you are presently reading this means that you have had !! knowledge of the CeCILL license and that you accept its terms. !! ========================================================================= #include SUBROUTINE PHOTO_INTERP( & zin ,& sin ,& vin ,& albin ,& t500in ,& t200in ,& ajout ) !---------------------------------------------------------------------- ! ... Loglinear interpolation for the photodissociation rates ! Note: this subroutine computes photorates for a vertical ! column at a given longitude and latitude ! This routine uses a six parameter table via a Taylor ! series expansion. ! Stacy Walters, Sep 30, 1996. Changed code to strictly limit ! the 200mb and 500mb temperature interpolation to the table ! endpoints; i.e. no extrapolation beyond the table is allowed. !---------------------------------------------------------------------- USE PHT_TABLES USE INCA_DIM IMPLICIT NONE !---------------------------------------------------------------------- ! ... Dummy arguments !---------------------------------------------------------------------- REAL, INTENT(in) :: zin(PLEV) ! pressure in millibars REAL, INTENT(in) :: sin ! secant solar zenith angle REAL, INTENT(in) :: vin(PLEV) ! o3 column density REAL, INTENT(in) :: albin(PLEV) ! surface albedo REAL, INTENT(in) :: t500in ! temp on 500mb surface REAL, INTENT(in) :: t200in ! temp on 200mb surface REAL, INTENT(out) :: ajout(jdim,PLEV) ! photodissociation rates !---------------------------------------------------------------------- ! ... Local variables !---------------------------------------------------------------------- INTEGER :: i, iz, is, iv, ial, nn, it500, it200 INTEGER :: i1, i2 INTEGER :: k INTEGER :: izl INTEGER :: index0 INTEGER :: base_ind INTEGER :: indv0(7,PLEV) INTEGER, DIMENSION(PLEV) :: altind, ratind, albind REAL :: v3std REAL :: psum REAL :: dels(6) REAL, DIMENSION(PLEV) :: v3rat, wght0, wght1, wght2, wght3, wght4 INTEGER :: first_level !---------------------------------------------------------------------- ! ... Find the zenith angle index ( same for all levels ) !---------------------------------------------------------------------- DO is = 1,zangdim IF( vsec(is) > sin ) THEN EXIT END IF END DO is = MAX( MIN( is,zangdim ) - 1,1 ) dels(2) = (sin - vsec(is)) * delang(is) !---------------------------------------------------------------------- ! ... Find the 500mb temp index ( same for all levels ) !---------------------------------------------------------------------- DO it500 = 1,t500dim IF( t500(it500) > t500in ) THEN EXIT END IF END DO it500 = MAX( MIN( it500,t500dim ) - 1,1 ) dels(5) = (t500in - t500(it500)) * delt500(it500) dels(5) = MAX( MIN( dels(5),1. ),0. ) !---------------------------------------------------------------------- ! ... Find the 200mb temp index ( same for all levels ) !---------------------------------------------------------------------- it200 = 1 dels(6) = (t200in - t200(it200)) * delt200(it200) dels(6) = MAX( MIN( dels(6),1. ),0. ) base_ind = is * offset(2) + it500 * offset(5) & + it200 * offset(6) - offset(7) izl = 1 DO k = PLEV,1,-1 !---------------------------------------------------------------------- ! ... Find albedo indicies !---------------------------------------------------------------------- DO ial = 1,albdim IF( albev(ial) > albin(k) ) THEN EXIT END IF END DO albind(k) = MAX( MIN( ial,albdim ) - 1,1 ) !---------------------------------------------------------------------- ! ... Find level indicies !---------------------------------------------------------------------- DO iz = izl,altdim IF( zz(iz) > zin(k) ) THEN izl = iz EXIT END IF END DO altind(k) = MAX( MIN( iz,altdim ) - 1,1 ) !---------------------------------------------------------------------- ! ... Find o3 ratio indicies !---------------------------------------------------------------------- i = MAX( MIN( altmxdim-1,INT( zin(k) ) ),0 ) v3std = vo3(i) + (zin(k) - FLOAT(i)) * (vo3(i+1) - vo3(i)) v3rat(k) = vin(k) / v3std DO iv = 1,o3ratdim IF( xv3(iv) > v3rat(k) ) THEN EXIT END IF END DO ratind(k) = MAX( MIN( iv,o3ratdim ) - 1,1 ) END DO !!! patch pour le 79 niveau - Fixed a bug in the last level (at 80km) !!! we put the last 3 level to the value of the fourth one !!! Anne et Didier 24 mai 2017 IF (PLEV .EQ. 79) THEN DO k = 1,3 iz = altind(4) iv = ratind(4) ial = albind(4) index0 = base_ind + ial * offset(4) & + iv * offset(3) + iz * offset(1) dels(1) = (zin(4) - zz(iz)) * delz(iz) dels(3) = (v3rat(4) - xv3(iv)) * delv(iv) dels(4) = (albin(4) - albev(ial)) * delalb(ial) !---------------------------------------------------------------------- ! ... Compute the weigths !---------------------------------------------------------------------- wght0(k) = 1. - SUM( dels ) wght1(k) = dels(1) wght3(k) = dels(3) wght4(k) = dels(4) !---------------------------------------------------------------------- ! ... Compute the indicies into ajl !---------------------------------------------------------------------- indv0(1,k) = index0 indv0(2,k) = index0 + offset(1) indv0(3,k) = index0 + offset(2) indv0(4,k) = index0 + offset(3) indv0(5,k) = index0 + offset(4) indv0(6,k) = index0 + offset(5) indv0(7,k) = index0 + offset(6) END DO first_level=4 ELSE first_level=1 ENDIF DO k = first_level,PLEV iz = altind(k) iv = ratind(k) ial = albind(k) index0 = base_ind + ial * offset(4) & + iv * offset(3) + iz * offset(1) dels(1) = (zin(k) - zz(iz)) * delz(iz) dels(3) = (v3rat(k) - xv3(iv)) * delv(iv) dels(4) = (albin(k) - albev(ial)) * delalb(ial) !---------------------------------------------------------------------- ! ... Compute the weigths !---------------------------------------------------------------------- wght0(k) = 1. - SUM( dels ) wght1(k) = dels(1) wght3(k) = dels(3) wght4(k) = dels(4) !---------------------------------------------------------------------- ! ... Compute the indicies into ajl !---------------------------------------------------------------------- indv0(1,k) = index0 indv0(2,k) = index0 + offset(1) indv0(3,k) = index0 + offset(2) indv0(4,k) = index0 + offset(3) indv0(5,k) = index0 + offset(4) indv0(6,k) = index0 + offset(5) indv0(7,k) = index0 + offset(6) END DO DO is = 1,jdim*PLEV k = MOD( is-1,PLEV ) + 1 nn = (is - 1)/PLEV + 1 psum = wght0(k) * ajl(indv0(1,k)+nn) & + wght1(k) * ajl(indv0(2,k)+nn) & + dels(2) * ajl(indv0(3,k)+nn) & + wght3(k) * ajl(indv0(4,k)+nn) & + wght4(k) * ajl(indv0(5,k)+nn) & + dels(5) * ajl(indv0(6,k)+nn) & + dels(6) * ajl(indv0(7,k)+nn) ajout(nn,k) = EXP( psum ) END DO END SUBROUTINE PHOTO_INTERP