!$Id: set_ub_vals.F90 157 2010-01-18 14:21:59Z 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: !! !! Didier Hauglustaine, LSCE, hauglustaine@cea.fr !! !! 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 #if defined(AERONLY) || defined(GES) SUBROUTINE XIOS_OXYDANT_READ(calday) !---------------------------------------------------------------------- ! ... Read OXYDANT distributions ! Didier Hauglustaine, IPSL, 1999. !---------------------------------------------------------------------- USE OXYDANT_COM USE MOD_GRID_INCA USE INCA_DIM USE MOD_INCA_PARA USE PRINT_INCA USE XIOS_INCA IMPLICIT NONE INCLUDE 'netcdf.inc' !---------------------------------------------------------------------- ! ... Dummy args !---------------------------------------------------------------------- REAL, INTENT(IN) :: calday !---------------------------------------------------------------------- ! ... Local variables !---------------------------------------------------------------------- INTEGER :: iret INTEGER :: varid INTEGER :: oldlotime, oldhitime INTEGER, DIMENSION(4) :: start3, count3 LOGICAL, SAVE :: first = .TRUE. !$OMP THREADPRIVATE(first) REAL, DIMENSION(12) :: timeo_inter REAL :: dt, dt1, tint INTEGER :: lonlev IF (first) THEN ALLOCATE(timeo(ntime_oxyd)) ALLOCATE(o3oxyd_year(PLON, PLEV, ntime_oxyd)) ALLOCATE(o1doxyd_year(PLON, PLEV, ntime_oxyd)) ALLOCATE(hno3oxyd_year(PLON, PLEV, ntime_oxyd)) ALLOCATE(ohoxyd_year(PLON, PLEV, ntime_oxyd)) ALLOCATE(h2o2oxyd_year(PLON, PLEV, ntime_oxyd)) ALLOCATE(no2oxyd_year(PLON, PLEV, ntime_oxyd)) ALLOCATE(no3oxyd_year(PLON, PLEV, ntime_oxyd)) ! write(lunout,*) '--------------' ! write(lunout,*) 'lecture oxyd : ' ! write(lunout,*) 'PLON = ', PLON ! write(lunout,*) 'PLEV = ', PLEV ! write(lunout,*) 'ntime_oxy = ', ntime_oxyd ! write(lunout,*) '--------------' CALL xios_inca_change_context("inca") CALL xios_inca_recv_field_glo("oxydtime_read", timeo ) CALL xios_inca_recv_field("O3_interp" , o3oxyd_year ) CALL xios_inca_recv_field("O1D_interp" , o1doxyd_year ) CALL xios_inca_recv_field("HNO3_interp" , hno3oxyd_year) CALL xios_inca_recv_field("OH_interp" , ohoxyd_year ) CALL xios_inca_recv_field("H2O2_interp" , h2o2oxyd_year) CALL xios_inca_recv_field("NO2_interp" , no2oxyd_year ) CALL xios_inca_recv_field("NO3_interp" , no3oxyd_year ) CALL xios_inca_send_field("NO3_oxyd_read", no3oxyd_year(:,:,12)) CALL xios_inca_change_context("LMDZ") first = .false. ENDIF lonlev = PLON * nlevo ! ... Check to see if model time still bounded by data set oldlotime = lotime oldhitime = hitime timeo_inter(:) = mod(timeo(:)/86400, 365.) CALL FINDPLB (timeo_inter, ntime_oxyd, calday, lotime) IF ( cyclical ) THEN hitime = MOD(lotime,ntime_oxyd) + 1 ELSE hitime = lotime + 1 END IF ! -------------------- interpolation sur le temps !---------------------------------------------------------------------- ! Note: 365 day year inconsistent with LMDz (360 days) !!! !---------------------------------------------------------------------- ! ... First interpolate linearly on current model time ! ! IF ( timeo(hitime) < timeo(lotime) ) THEN dt = 365. - mod(timeo(lotime)/(86400),365.) + mod(timeo(hitime)/(86400),365.) IF (calday <= mod(timeo(hitime)/(86400),365.) ) THEN dt1 = 365. - mod(timeo(lotime)/(86400),365.) + calday ELSE dt1 = calday - mod(timeo(lotime)/(86400),365.) END IF ELSE dt = mod(timeo(hitime)/(86400),365.) - mod(timeo(lotime)/(86400),365.) dt1 = calday - mod(timeo(lotime)/(86400),365.) END IF tint = dt1/dt ohoxyd(:,:) = ohoxyd_year(:,:,lotime) + (ohoxyd_year(:,:,hitime)-ohoxyd_year(:,:,lotime))*tint o1doxyd(:,:) = o1doxyd_year(:,:,lotime) + (o1doxyd_year(:,:,hitime)-o1doxyd_year(:,:,lotime))*tint o3oxyd(:,:) = o3oxyd_year(:,:,lotime) + (o3oxyd_year(:,:,hitime)-o3oxyd_year(:,:,lotime))*tint no3oxyd(:,:) = no3oxyd_year(:,:,lotime) + (no3oxyd_year(:,:,hitime)-no3oxyd_year(:,:,lotime))*tint h2o2oxyd(:,:) = h2o2oxyd_year(:,:,lotime) + (h2o2oxyd_year(:,:,hitime)-h2o2oxyd_year(:,:,lotime))*tint hno3oxyd(:,:) = hno3oxyd_year(:,:,lotime) + (hno3oxyd_year(:,:,hitime)-hno3oxyd_year(:,:,lotime))*tint no2oxyd(:,:) = no2oxyd_year(:,:,lotime) + (no2oxyd_year(:,:,hitime)-no2oxyd_year(:,:,lotime))*tint CALL xios_inca_change_context("inca") CALL xios_inca_send_field("O3_oxyd" , o3oxyd ) CALL xios_inca_send_field("O1D_oxyd" , o1doxyd ) CALL xios_inca_send_field("HNO3_oxyd" , hno3oxyd ) CALL xios_inca_send_field("OH_oxyd" , ohoxyd ) CALL xios_inca_send_field("H2O2_oxyd" , h2o2oxyd ) CALL xios_inca_send_field("NO2_oxyd" , no2oxyd ) CALL xios_inca_send_field("NO3_oxyd" , no3oxyd ) CALL xios_inca_change_context("LMDZ") END SUBROUTINE XIOS_OXYDANT_READ #else ! IF not defined aeronly and ges SUBROUTINE OZCLIM_READ(calday) !---------------------------------------------------------------------- ! ... Read ozone distributions ! Didier Hauglustaine, IPSL, 1999. !---------------------------------------------------------------------- USE O3CLIM_COM USE INCA_DIM USE MOD_INCA_PARA USE PRINT_INCA IMPLICIT NONE INCLUDE 'netcdf.inc' !---------------------------------------------------------------------- ! ... Dummy args !---------------------------------------------------------------------- REAL, INTENT(IN) :: calday !---------------------------------------------------------------------- ! ... Local variables !---------------------------------------------------------------------- INTEGER :: iret INTEGER :: varid INTEGER :: oldlotime, oldhitime INTEGER, DIMENSION(3) :: start, count REAL :: o3climbd_glo(nbp_glo,nlevc,2) ! ... Check to see if model time still bounded by data set oldlotime = lotime oldhitime = hitime CALL FINDPLB (timec, ntimec, calday, lotime) IF ( cyclical ) THEN hitime = MOD(lotime,ntimec) + 1 ELSE hitime = lotime + 1 END IF ! ... Read new data if needed IF ( hitime /= oldhitime ) THEN loind = hiind hiind = MOD( loind, 2) + 1 !$OMP MASTER IF (is_mpi_root) THEN iret = nf_inq_varid(ncidc, 'O3', varid) CALL check_err(iret, 'OZCLIM_READ','problem to check varid O3') count(1) = nbp_glo count(2) = nlevc count(3) = 1 start(1) = 1 start(2) = 1 WRITE(lunout,*) 'OZCLIM_READ : read new data for time ', timec(hitime) start(3) = hitime iret = nf_get_vara_double(ncidc, varid, start, count,o3climbd_glo(1,1,hiind)) CALL check_err(iret, 'OZCLIM_READ', 'problem to read O3 hiind') ENDIF !$OMP END MASTER CALL scatter(o3climbd_glo(:,:,hiind),o3climbd(:,:,hiind)) IF ( lotime /= oldhitime ) THEN !$OMP MASTER IF (is_mpi_root) THEN WRITE(lunout,*) 'OZCLIM_READ : read new data for time ', timec(lotime) start(3) = lotime iret = nf_get_vara_double(ncidc, varid, start, count, o3climbd_glo(1,1,loind)) CALL check_err(iret, 'OZCLIM_READ', 'problem to read O3 loind') ENDIF !$OMP END MASTER CALL scatter(o3climbd_glo(:,:,loind),o3climbd(:,:,loind)) END IF END IF END SUBROUTINE OZCLIM_READ SUBROUTINE OZCLIM_INTERP(calday, pmid) !---------------------------------------------------------------------- ! ... Interpolate ozone on model time and on pressure levels ! Didier Hauglustaine, IPSL, 1999. !---------------------------------------------------------------------- USE O3CLIM_COM USE INCA_DIM IMPLICIT NONE !---------------------------------------------------------------------- ! ... Dummy args !---------------------------------------------------------------------- REAL, INTENT(IN) :: calday REAL, INTENT(IN) :: pmid(PLON,PLEV) !---------------------------------------------------------------------- ! ... Local variables !---------------------------------------------------------------------- REAL :: dt, dt1, tint REAL :: rma = 48./28. REAL :: ploc(klonc,nlevc) INTEGER :: lonlev INTEGER :: i INTEGER :: index(PLON,PLEV) index = -9999 lonlev = klonc * nlevc DO i = 1, klonc ploc(i,:) = presc(:) END DO !---------------------------------------------------------------------- ! Note: 365 day year inconsistent with LMDz (360 days) !!! !---------------------------------------------------------------------- ! ... First interpolate linearly on current model time IF ( timec(hitime) < timec(lotime) ) THEN dt = 365. - timec(lotime) + timec(hitime) IF (calday <= timec(hitime) ) THEN dt1 = 365. - timec(lotime) + calday ELSE dt1 = calday - timec(lotime) END IF ELSE dt = timec(hitime) - timec(lotime) dt1 = calday - timec(lotime) END IF tint = dt1/dt CALL linintp (lonlev, & 0., & 1., & tint, & o3climbd(1,1,loind), & o3climbd(1,1,hiind), & o3climcr) ! ... Put on model pressure levels CALL pinterp (PLON, & o3climcr, & ploc, & nlevc, & o3clim, & pmid, & PLEV, & index, & .false. ) ! ... vmr to mmr o3clim = o3clim * rma END SUBROUTINE OZCLIM_INTERP SUBROUTINE OZRESET (mmr, pmid, temp) !---------------------------------------------------------------------- ! ... Reset O3 boundary conditions for use in radiation ! Didier Hauglustaine, IPSL, 2002. !---------------------------------------------------------------------- USE O3CLIM_COM USE SPECIES_NAMES USE OXYDANT_COM USE AIRPLANE_SRC, ONLY : itrop USE CHEM_MODS, ONLY : adv_mass USE INCA_DIM IMPLICIT NONE !---------------------------------------------------------------------- ! ... Dummy args !---------------------------------------------------------------------- REAL, INTENT(IN) :: pmid(PLON,PLEV) REAL, INTENT(IN) :: temp(PLON,PLEV) REAL, INTENT(INOUT) :: mmr(PLON,PLEV,PCNST) !---------------------------------------------------------------------- ! ... Local variables !---------------------------------------------------------------------- INTEGER :: i, k INTEGER :: ktrop, krelax, k380 REAL, PARAMETER :: dry_mass = 28.966 REAL :: theta(PLON,PLEV) theta(:,:) = temp(:,:)*(1.e5/pmid(:,:))**0.286 DO i = 1, PLON ktrop = MAX(1,itrop(i)) DO k = ktrop, PLEV IF (theta(i,k) >= 380.) THEN k380 = k EXIT ENDIF ENDDO ! krelax = MIN(ktrop+2,PLEV) ! 2 levels above tropopause krelax = k380 ! 380K theta surface mmr(i,krelax:PLEV,id_O3) = o3clim(i,krelax:PLEV) END DO END SUBROUTINE OZRESET SUBROUTINE FIELD_PREP (mmr) !---------------------------------------------------------------------- ! ... Prepare fields to be used in LMDz radiation ! Didier Hauglustaine, IPSL, 2002. 2018. !---------------------------------------------------------------------- USE SPECIES_NAMES USE CHEM_MODS, ONLY : o3_inca, h2o_inca, ch4_inca, n2o_inca, cfc11_inca, cfc12_inca, & ch4_inca_2d, n2o_inca_2d, cfc11_inca_2d, cfc12_inca_2d, adv_mass USE INCA_DIM IMPLICIT NONE !---------------------------------------------------------------------- ! ... Dummy args !---------------------------------------------------------------------- REAL, INTENT(IN) :: mmr(PLON,PLEV,PCNST) !---------------------------------------------------------------------- ! ... Local variables !---------------------------------------------------------------------- INTEGER :: jj, jk, jki, jkf INTEGER, PARAMETER :: ng1 = 2, ng1p1 = ng1 + 1 !From raddimlw.h ! Fields to be used in LMDz in mass mixing ratio #ifdef NMHC o3_inca(:,:) = mmr(:,:,id_O3) #endif ! ch4_inca(:,:) = mmr(:,:,id_CH4) ! n2o_inca(:,:) = mmr(:,:,id_N2O) ! cfc11_inca(:,:) = mmr(:,:,id_CFC11) ! cfc12_inca(:,:) = mmr(:,:,id_CFC12) ! h2o_inca(:,:) = mmr(:,:,id_H2O) ! ! ! ! ! Ozone field POZON (kg/kg) ! ! Note: For POZON (used in radlwsw) simply use o3_inca from above (in mmr) ! ! ! 2D fields for CH4, N2O, CFC11 and CFC12 (in mmr) ! DO jj = 1 , PLEV ! jki=(jj-1)*ng1p1+1 ! jkf=jki+ng1 ! DO jk = jki, jkf ! ch4_inca_2d(:,jk) = ch4_inca(:,jj) ! n2o_inca_2d(:,jk) = n2o_inca(:,jj) ! cfc11_inca_2d(:,jk) = cfc11_inca(:,jj) ! cfc12_inca_2d(:,jk) = cfc12_inca(:,jj) ! ENDDO ! ENDDO END SUBROUTINE FIELD_PREP SUBROUTINE OZLIN_READ(calday) !---------------------------------------------------------------------- ! ... Read ozone linear coefficients ! Didier Hauglustaine, IPSL, 1999. !---------------------------------------------------------------------- USE O3LIN_COM USE INCA_DIM USE MOD_INCA_PARA USE PRINT_INCA IMPLICIT NONE INCLUDE 'netcdf.inc' !---------------------------------------------------------------------- ! ... Dummy args !---------------------------------------------------------------------- REAL, INTENT(IN) :: calday !---------------------------------------------------------------------- ! ... Local variables !---------------------------------------------------------------------- INTEGER :: iret INTEGER :: varid1,varid2,varid3,varid4,varid5 INTEGER :: varid6,varid7,varid8,varid9 INTEGER :: oldlotime, oldhitime INTEGER, DIMENSION(3) :: start, count REAL :: A1bd_glo(nbp_glo,nlevl,2) REAL :: A2bd_glo(nbp_glo,nlevl,2) REAL :: A3bd_glo(nbp_glo,nlevl,2) REAL :: A4bd_glo(nbp_glo,nlevl,2) REAL :: A5bd_glo(nbp_glo,nlevl,2) REAL :: A6bd_glo(nbp_glo,nlevl,2) REAL :: A7bd_glo(nbp_glo,nlevl,2) REAL :: A8bd_glo(nbp_glo,nlevl,2) REAL :: A9bd_glo(nbp_glo,nlevl,2) ! ... Check to see if model time still bounded by data set oldlotime = lotime oldhitime = hitime CALL FINDPLB (timel, ntimel, calday, lotime) IF ( cyclical ) THEN hitime = MOD(lotime,ntimel) + 1 ELSE hitime = lotime + 1 END IF ! ... Read new data if needed IF ( hitime /= oldhitime ) THEN loind = hiind hiind = MOD( loind, 2) + 1 !$OMP MASTER IF (is_mpi_root) THEN iret = nf_inq_varid(ncidl, 'A1', varid1) CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A1') iret = nf_inq_varid(ncidl, 'A2', varid2) CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A2') iret = nf_inq_varid(ncidl, 'A3', varid3) CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A3') iret = nf_inq_varid(ncidl, 'A4', varid4) CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A4') iret = nf_inq_varid(ncidl, 'A5', varid5) CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A5') iret = nf_inq_varid(ncidl, 'A6', varid6) CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A6') iret = nf_inq_varid(ncidl, 'A7', varid7) CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A7') iret = nf_inq_varid(ncidl, 'A8', varid8) CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A8') iret = nf_inq_varid(ncidl, 'A9', varid9) CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A9') count(1) = nbp_glo count(2) = nlevl count(3) = 1 start(1) = 1 start(2) = 1 WRITE(lunout,*) 'OZLIN_READ : read new data for time ', timel(hitime) start(3) = hitime iret = nf_get_vara_double(ncidl, varid1, start, count, A1bd_glo(1,1,hiind)) CALL check_err(iret, 'OCLIN_READ', 'problem to read A1bd for hiind') iret = nf_get_vara_double(ncidl, varid2, start, count, A2bd_glo(1,1,hiind)) CALL check_err(iret, 'OCLIN_READ', 'problem to read A2bd for hiind') iret = nf_get_vara_double(ncidl, varid3, start, count, A3bd_glo(1,1,hiind)) CALL check_err(iret, 'OCLIN_READ', 'problem to read A3bd for hiind') iret = nf_get_vara_double(ncidl, varid4, start, count, A4bd_glo(1,1,hiind)) CALL check_err(iret, 'OCLIN_READ', 'problem to read A4bd for hiind') iret = nf_get_vara_double(ncidl, varid5, start, count, A5bd_glo(1,1,hiind)) CALL check_err(iret, 'OCLIN_READ', 'problem to read A5bd for hiind') iret = nf_get_vara_double(ncidl, varid6, start, count, A6bd_glo(1,1,hiind)) CALL check_err(iret, 'OCLIN_READ', 'problem to read A6bd for hiind') iret = nf_get_vara_double(ncidl, varid7, start, count, A7bd_glo(1,1,hiind)) CALL check_err(iret, 'OCLIN_READ', 'problem to read A7bd for hiind') iret = nf_get_vara_double(ncidl, varid8, start, count, A8bd_glo(1,1,hiind)) CALL check_err(iret, 'OCLIN_READ', 'problem to read A8bd for hiind') iret = nf_get_vara_double(ncidl, varid9, start, count, A9bd_glo(1,1,hiind)) CALL check_err(iret, 'OCLIN_READ', 'problem to read A9bd for hiind') ENDIF !$OMP END MASTER CALL scatter(A1bd_glo(:,:,hiind),A1bd(:,:,hiind)) CALL scatter(A2bd_glo(:,:,hiind),A2bd(:,:,hiind)) CALL scatter(A3bd_glo(:,:,hiind),A3bd(:,:,hiind)) CALL scatter(A4bd_glo(:,:,hiind),A4bd(:,:,hiind)) CALL scatter(A5bd_glo(:,:,hiind),A5bd(:,:,hiind)) CALL scatter(A6bd_glo(:,:,hiind),A6bd(:,:,hiind)) CALL scatter(A7bd_glo(:,:,hiind),A7bd(:,:,hiind)) CALL scatter(A8bd_glo(:,:,hiind),A8bd(:,:,hiind)) CALL scatter(A9bd_glo(:,:,hiind),A9bd(:,:,hiind)) IF ( lotime /= oldhitime ) THEN !$OMP MASTER IF (is_mpi_root) THEN WRITE(lunout,*) 'OZLIN_READ : read new data for time ', timel(lotime) start(3) = lotime iret = nf_get_vara_double(ncidl, varid1, start, count, A1bd_glo(1,1,loind)) CALL check_err(iret, 'OZLIN_READ', 'problem to read A1bd for loind') iret = nf_get_vara_double(ncidl, varid2, start, count, A2bd_glo(1,1,loind)) CALL check_err(iret, 'OZLIN_READ', 'problem to read A2bd for loind') iret = nf_get_vara_double(ncidl, varid3, start, count, A3bd_glo(1,1,loind)) CALL check_err(iret, 'OZLIN_READ', 'problem to read A3bd for loind') iret = nf_get_vara_double(ncidl, varid4, start, count, A4bd_glo(1,1,loind)) CALL check_err(iret, 'OZLIN_READ', 'problem to read A4bd for loind') iret = nf_get_vara_double(ncidl, varid5, start, count, A5bd_glo(1,1,loind)) CALL check_err(iret, 'OZLIN_READ', 'problem to read A5bd for loind') iret = nf_get_vara_double(ncidl, varid6, start, count, A6bd_glo(1,1,loind)) CALL check_err(iret, 'OZLIN_READ', 'problem to read A6bd for loind') iret = nf_get_vara_double(ncidl, varid7, start, count, A7bd_glo(1,1,loind)) CALL check_err(iret, 'OZLIN_READ', 'problem to read A7bd for loind') iret = nf_get_vara_double(ncidl, varid8, start, count, A8bd_glo(1,1,loind)) CALL check_err(iret, 'OZLIN_READ', 'problem to read A8bd for loind') iret = nf_get_vara_double(ncidl, varid9, start, count, A9bd_glo(1,1,loind)) CALL check_err(iret, 'OZLIN_READ', 'problem to read A9bd for loind') ENDIF !$OMP END MASTER CALL scatter(A1bd_glo(:,:,loind),A1bd(:,:,loind)) CALL scatter(A2bd_glo(:,:,loind),A2bd(:,:,loind)) CALL scatter(A3bd_glo(:,:,loind),A3bd(:,:,loind)) CALL scatter(A4bd_glo(:,:,loind),A4bd(:,:,loind)) CALL scatter(A5bd_glo(:,:,loind),A5bd(:,:,loind)) CALL scatter(A6bd_glo(:,:,loind),A6bd(:,:,loind)) CALL scatter(A7bd_glo(:,:,loind),A7bd(:,:,loind)) CALL scatter(A8bd_glo(:,:,loind),A8bd(:,:,loind)) CALL scatter(A9bd_glo(:,:,loind),A9bd(:,:,loind)) END IF END IF END SUBROUTINE OZLIN_READ SUBROUTINE OZLIN_INTERP(calday, pmid) !---------------------------------------------------------------------- ! ... Interpolate ozone linear coefficients on model time and on pressure levels !R. Valorso !---------------------------------------------------------------------- USE O3LIN_COM USE INCA_DIM IMPLICIT NONE !---------------------------------------------------------------------- ! ... Dummy args !---------------------------------------------------------------------- REAL, INTENT(IN) :: calday REAL, INTENT(IN) :: pmid(PLON,PLEV) !---------------------------------------------------------------------- ! ... Local variables !---------------------------------------------------------------------- REAL :: dt, dt1, tint REAL :: ploc(klonl,nlevl) INTEGER :: lonlev INTEGER :: i INTEGER :: index(PLON,PLEV) index = -9999 lonlev = klonl * nlevl DO i = 1, klonl ploc(i,:) = presl(:) END DO !---------------------------------------------------------------------- ! Note: 365 day year inconsistent with LMDz (360 days) !!! !---------------------------------------------------------------------- ! ... First interpolate linearly on current model time IF ( timel(hitime) < timel(lotime) ) THEN dt = 365. - timel(lotime) + timel(hitime) IF (calday <= timel(hitime) ) THEN dt1 = 365. - timel(lotime) + calday ELSE dt1 = calday - timel(lotime) END IF ELSE dt = timel(hitime) - timel(lotime) dt1 = calday - timel(lotime) END IF tint = dt1/dt CALL linintp (lonlev, & 0., & 1., & tint, & A1bd(1,1,loind), & A1bd(1,1,hiind), & A1cr) ! ... Put on model pressure levels CALL pinterp (PLON, & A1cr, & ploc, & nlevl, & A1, & pmid, & PLEV, & index, & .false. ) ! ... CALL linintp (lonlev, & 0., & 1., & tint, & A2bd(1,1,loind), & A2bd(1,1,hiind), & A2cr) ! ... Put on model pressure levels index = -9999 CALL pinterp (PLON, & A2cr, & ploc, & nlevl, & A2, & pmid, & PLEV, & index, & .false. ) ! ... CALL linintp (lonlev, & 0., & 1., & tint, & A3bd(1,1,loind), & A3bd(1,1,hiind), & A3cr) ! ... Put on model pressure levels index = -9999 CALL pinterp (PLON, & A3cr, & ploc, & nlevl, & A3, & pmid, & PLEV, & index, & .false. ) ! ... CALL linintp (lonlev, & 0., & 1., & tint, & A4bd(1,1,loind), & A4bd(1,1,hiind), & A4cr) ! ... Put on model pressure levels index = -9999 CALL pinterp (PLON, & A4cr, & ploc, & nlevl, & A4, & pmid, & PLEV, & index, & .false. ) ! ... CALL linintp (lonlev, & 0., & 1., & tint, & A5bd(1,1,loind), & A5bd(1,1,hiind), & A5cr) ! ... Put on model pressure levels index = -9999 CALL pinterp (PLON, & A5cr, & ploc, & nlevl, & A5, & pmid, & PLEV, & index, & .false. ) ! ... CALL linintp (lonlev, & 0., & 1., & tint, & A6bd(1,1,loind), & A6bd(1,1,hiind), & A6cr) ! ... Put on model pressure levels index = -9999 CALL pinterp (PLON, & A6cr, & ploc, & nlevl, & A6, & pmid, & PLEV, & index, & .false. ) ! ... CALL linintp (lonlev, & 0., & 1., & tint, & A7bd(1,1,loind), & A7bd(1,1,hiind), & A7cr) ! ... Put on model pressure levels index = -9999 CALL pinterp (PLON, & A7cr, & ploc, & nlevl, & A7, & pmid, & PLEV, & index, & .false. ) ! ... CALL linintp (lonlev, & 0., & 1., & tint, & A8bd(1,1,loind), & A8bd(1,1,hiind), & A8cr) ! ... Put on model pressure levels index = -9999 CALL pinterp (PLON, & A8cr, & ploc, & nlevl, & A8, & pmid, & PLEV, & index, & .false. ) ! ... CALL linintp (lonlev, & 0., & 1., & tint, & A9bd(1,1,loind), & A9bd(1,1,hiind), & A9cr) ! ... Put on model pressure levels index = -9999 CALL pinterp (PLON, & A9cr, & ploc, & nlevl, & A9, & pmid, & PLEV, & index, & .false. ) END SUBROUTINE OZLIN_INTERP #ifdef STRAT SUBROUTINE SAD_READ(calday) !---------------------------------------------------------------------- ! ... Read sulfate data from CMIP5 ! R. Valorso !---------------------------------------------------------------------- USE SAD_COM USE INCA_DIM USE MOD_INCA_PARA USE PRINT_INCA IMPLICIT NONE INCLUDE 'netcdf.inc' !---------------------------------------------------------------------- ! ... Dummy args !---------------------------------------------------------------------- REAL, INTENT(IN) :: calday !---------------------------------------------------------------------- ! ... Local variables !---------------------------------------------------------------------- INTEGER :: iret INTEGER :: varid1,varid2,varid3,varid4,varid5 INTEGER :: oldlotime, oldhitime INTEGER, DIMENSION(3) :: start, count REAL :: PRES_bd_glo(nbp_glo,nlevs,2) REAL :: SAD_bd_glo(nbp_glo,nlevs,2) REAL :: MASS_bd_glo(nbp_glo,nlevs,2) REAL :: VOL_bd_glo(nbp_glo,nlevs,2) REAL :: RMEAN_bd_glo(nbp_glo,nlevs,2) ! ... Check to see if model time still bounded by data set oldlotime = lotime oldhitime = hitime CALL FINDPLB (times, ntimes, calday, lotime) IF ( cyclical ) THEN hitime = MOD(lotime,ntimes) + 1 ELSE hitime = lotime + 1 END IF ! ... Read new data if needed IF ( hitime /= oldhitime ) THEN loind = hiind hiind = MOD( loind, 2) + 1 !$OMP MASTER IF (is_mpi_root) THEN iret = nf_inq_varid(ncids, 'PRES', varid1) CALL check_err(iret, 'SAF_READ', 'problem to get id for pres') iret = nf_inq_varid(ncids, 'SAD', varid2) CALL check_err(iret, 'SAD_READ', 'problem to get id for sad') iret = nf_inq_varid(ncids, 'MASS', varid3) CALL check_err(iret, 'SAD_READ', 'problem to get id for mass') iret = nf_inq_varid(ncids, 'VOL', varid4) CALL check_err(iret, 'SAD_READ', 'problem to get id for vol') iret = nf_inq_varid(ncids, 'RMEAN', varid5) CALL check_err(iret, 'SAD_READ', 'problem to get id for rmean') count(1) = nbp_glo count(2) = nlevs count(3) = 1 start(1) = 1 start(2) = 1 start(3) = hitime iret = nf_get_vara_double(ncids, varid1, start, count, PRES_bd_glo(1,1,hiind)) CALL check_err(iret, 'SAD_READ', 'problem to read PRES_bd hiind') iret = nf_get_vara_double(ncids, varid2, start, count, SAD_bd_glo(1,1,hiind)) CALL check_err(iret, 'SAD_READ', 'problem to read SAD_bd hiind') iret = nf_get_vara_double(ncids, varid3, start, count, MASS_bd_glo(1,1,hiind)) CALL check_err(iret, 'SAD_READ', 'problem to read MASS_bd hiind') iret = nf_get_vara_double(ncids, varid4, start, count, VOL_bd_glo(1,1,hiind)) CALL check_err(iret, 'SAD_READ', 'problem to read VOL_bd hiind') iret = nf_get_vara_double(ncids, varid5, start, count, RMEAN_bd_glo(1,1,hiind)) CALL check_err(iret, 'SAD_READ', 'problem to read RMEAN_bd hiind') ENDIF !$OMP END MASTER CALL scatter(PRES_bd_glo(:,:,hiind),PRES_bd(:,:,hiind)) CALL scatter(SAD_bd_glo(:,:,hiind),SAD_bd(:,:,hiind)) CALL scatter(MASS_bd_glo(:,:,hiind),MASS_bd(:,:,hiind)) CALL scatter(VOL_bd_glo(:,:,hiind),VOL_bd(:,:,hiind)) CALL scatter(RMEAN_bd_glo(:,:,hiind),RMEAN_bd(:,:,hiind)) IF ( lotime /= oldhitime ) THEN !$OMP MASTER IF (is_mpi_root) THEN start(3) = lotime iret = nf_get_vara_double(ncids, varid1, start, count, PRES_bd_glo(1,1,loind)) CALL check_err(iret, 'SAD_READ', 'problem to read PRES_bd loind') iret = nf_get_vara_double(ncids, varid2, start, count, SAD_bd_glo(1,1,loind)) CALL check_err(iret, 'SAD_READ', 'problem to read SAD_bd loind') iret = nf_get_vara_double(ncids, varid3, start, count, MASS_bd_glo(1,1,loind)) CALL check_err(iret, 'SAD_READ', 'problem to read MASS_bd loind') iret = nf_get_vara_double(ncids, varid4, start, count, VOL_bd_glo(1,1,loind)) CALL check_err(iret, 'SAD_READ', 'problem to read VOL_bd loind') iret = nf_get_vara_double(ncids, varid5, start, count, RMEAN_bd_glo(1,1,loind)) CALL check_err(iret, 'SAD_READ', 'problem to read RMEAN_bd loind') ENDIF !$OMP END MASTER CALL scatter(PRES_bd_glo(:,:,loind),PRES_bd(:,:,loind)) CALL scatter(SAD_bd_glo(:,:,loind),SAD_bd(:,:,loind)) CALL scatter(MASS_bd_glo(:,:,loind),MASS_bd(:,:,loind)) CALL scatter(VOL_bd_glo(:,:,loind),VOL_bd(:,:,loind)) CALL scatter(RMEAN_bd_glo(:,:,loind),RMEAN_bd(:,:,loind)) END IF END IF END SUBROUTINE SAD_READ SUBROUTINE SAD_INTERP(calday, pmid) !---------------------------------------------------------------------- ! ... Interpolate sulfate data from CMIP 5 on model time and on pressure levels ! R. Valorso !---------------------------------------------------------------------- USE SAD_COM USE INCA_DIM IMPLICIT NONE !---------------------------------------------------------------------- ! ... Dummy args !---------------------------------------------------------------------- REAL, INTENT(IN) :: calday REAL, INTENT(IN) :: pmid(PLON,PLEV) !---------------------------------------------------------------------- ! ... Local variables !---------------------------------------------------------------------- REAL :: dt, dt1, tint REAL :: ploc(klons,nlevs) INTEGER :: lonlev INTEGER :: i INTEGER :: index(PLON,PLEV) index = -9999 lonlev = klons * nlevs !---------------------------------------------------------------------- ! ... First interpolate linearly on current model time IF ( times(hitime) < times(lotime) ) THEN dt = 365. - times(lotime) + times(hitime) IF (calday <= times(hitime) ) THEN dt1 = 365. - times(lotime) + calday ELSE dt1 = calday - times(lotime) END IF ELSE dt = times(hitime) - times(lotime) dt1 = calday - times(lotime) END IF tint = dt1/dt CALL linintp (lonlev, & 0., & 1., & tint, & PRES_bd(1,1,loind), & PRES_bd(1,1,hiind), & PRES_cr) DO i = 1, klons ploc(i,:) = PRES_cr(i,:) ENDDO ! ... CALL linintp (lonlev, & 0., & 1., & tint, & SAD_bd(1,1,loind), & SAD_bd(1,1,hiind), & SAD_cr) ! ... Put on model pressure levels index = -9999 CALL pinterp (PLON, & SAD_cr, & ploc, & nlevs, & SAD, & pmid, & PLEV, & index, & .false. ) ! ... CALL linintp (lonlev, & 0., & 1., & tint, & MASS_bd(1,1,loind), & MASS_bd(1,1,hiind), & MASS_cr) ! ... Put on model pressure levels index = -9999 CALL pinterp (PLON, & MASS_cr, & ploc, & nlevs, & MASS, & pmid, & PLEV, & index, & .false. ) ! ... CALL linintp (lonlev, & 0., & 1., & tint, & VOL_bd(1,1,loind), & VOL_bd(1,1,hiind), & VOL_cr) ! ... Put on model pressure levels index = -9999 CALL pinterp (PLON, & VOL_cr, & ploc, & nlevs, & VOL, & pmid, & PLEV, & index, & .false. ) ! ... CALL linintp (lonlev, & 0., & 1., & tint, & RMEAN_bd(1,1,loind), & RMEAN_bd(1,1,hiind), & RMEAN_cr) ! ... Put on model pressure levels index = -9999 CALL pinterp (PLON, & RMEAN_cr, & ploc, & nlevs, & RMEAN, & pmid, & PLEV, & index, & .false. ) END SUBROUTINE SAD_INTERP #endif #endif #if !defined(GES) && !defined(DUSS) SUBROUTINE BOUNDSPC (dtinv, mmr, pmid, temp) !---------------------------------------------------------------------- ! ... Impose boundary conditions for chemical tracers ! Didier Hauglustaine, IPSL, 1999. !---------------------------------------------------------------------- USE O3CLIM_COM USE SPECIES_NAMES USE OXYDANT_COM USE AIRPLANE_SRC, ONLY : itrop USE CHEM_MODS, ONLY : adv_mass, invariants USE INCA_DIM #ifdef STRAT USE LGLIVED_MOD #endif USE PARAM_CHEM USE RATE_INDEX_MOD IMPLICIT NONE !---------------------------------------------------------------------- ! ... Dummy args !---------------------------------------------------------------------- REAL, INTENT(IN) :: dtinv REAL, INTENT(IN) :: pmid(PLON,PLEV) REAL, INTENT(IN) :: temp(PLON,PLEV) REAL, INTENT(INOUT) :: mmr(PLON,PLEV,PCNST) !---------------------------------------------------------------------- ! ... Local variables !---------------------------------------------------------------------- INTEGER :: i, k INTEGER :: ktrop, krelax, k380 REAL :: taurelax, facrelax REAL :: delt REAL, PARAMETER :: dry_mass = 28.966 REAL :: theta(PLON,PLEV) delt = 1./dtinv taurelax = 86400.*10. ! (10 days) facrelax = 1. - EXP( -delt / taurelax ) theta(:,:) = temp(:,:)*(1.e5/pmid(:,:))**0.286 DO i = 1, PLON ktrop = itrop(i) DO k = ktrop, PLEV IF (theta(i,k) >= 380.) THEN k380 = k EXIT ENDIF ENDDO ! krelax = MIN(ktrop+2,PLEV) ! 2 levels above tropopause krelax = k380 ! 380K theta surface ! ... Impose NOy to climatologies ! mmr(i,krelax:PLEV,id_NO) = nomzrt(i,krelax:PLEV) * adv_mass(id_NO)/dry_mass #ifdef AERONLY mmr(i,krelax:PLEV,id_NO2) = invariants(i,krelax:PLEV, inv_NO2INV)/invariants(i,krelax:PLEV, inv_M) * adv_mass(id_NO2)/dry_mass mmr(i,krelax:PLEV,id_HNO3) = invariants(i,krelax:PLEV, inv_HNO3INV)/invariants(i,krelax:PLEV, inv_M) * adv_mass(id_HNO3)/dry_mass #else ! ... Relax O3 to climatologies if (trim(flag_o3) .eq. "o3clim") then mmr(i,krelax:PLEV,id_O3) = o3clim(i,krelax:PLEV) elseif (trim(flag_o3) .eq. "o3lin") then mmr(i,krelax:PLEV,id_O3) = mmr(i,ktrop:PLEV,id_O3L) endif ! ... Impose inert O3 to O3 above tropopause mmr(i,ktrop:PLEV,id_O3I) = mmr(i,ktrop:PLEV,id_O3) ! ... Impose stratospheric O3 to O3 above tropopause mmr(i,ktrop:PLEV,id_O3S) = mmr(i,ktrop:PLEV,id_O3) #endif ! ... Impose CH4 ! mmr(i,:,id_CH4) = 1760E-9 * 16./dry_mass ! forcage CH4 pour IPCC2000 ! mmr(i,:,id_CH4) = 791.6E-9 * 16./dry_mass ! forcage CH4 pour IPCC2000 #ifdef STRAT ! Stratospheric Version ! Impose long_lived species at the surface CMIP5 mmr(i,1,id_CH4) = conc_cfc(1) * adv_mass(id_CH4) / dry_mass mmr(i,1,id_N2O) = conc_cfc(2) * adv_mass(id_N2O) / dry_mass mmr(i,1,id_CFC11) = conc_cfc(3) * adv_mass(id_CFC11) / dry_mass mmr(i,1,id_CFC12) = conc_cfc(4) * adv_mass(id_CFC12) / dry_mass mmr(i,1,id_CFC113) = conc_cfc(5) * adv_mass(id_CFC113) / dry_mass mmr(i,1,id_HCFC22) = conc_cfc(6) * adv_mass(id_HCFC22) / dry_mass mmr(i,1,id_HFC134a) = conc_cfc(7) * adv_mass(id_HFC134a) / dry_mass mmr(i,1,id_HCFC141) = conc_cfc(8) * adv_mass(id_HCFC141) / dry_mass mmr(i,1,id_HCFC142) = conc_cfc(9) * adv_mass(id_HCFC142) / dry_mass mmr(i,1,id_Ha1301) = conc_cfc(10) * adv_mass(id_Ha1301) / dry_mass mmr(i,1,id_Ha1211) = conc_cfc(11) * adv_mass(id_Ha1211) / dry_mass mmr(i,1,id_MCF) = conc_cfc(12) * adv_mass(id_MCF) / dry_mass mmr(i,1,id_CCl4) = conc_cfc(13) * adv_mass(id_CCl4) / dry_mass mmr(i,1,id_CH3Cl) = conc_cfc(14) * adv_mass(id_CH3Cl) / dry_mass mmr(i,1,id_CH3Br) = conc_cfc(15) * adv_mass(id_CH3Br) / dry_mass #endif END DO END SUBROUTINE BOUNDSPC #endif SUBROUTINE FDTROPOPAUSE (nx, ny, nz, pmid, temp) !**************************************************************************** ! Purpose ! ------- ! ! Calculates 2D fields of thermal tropopause pressures and tropopause ! level indices for given 3D fields of temperature and pressure. ! ! Methods ! ------- ! - Subroutine stattrop calculates the thermal tropopause pressure (Pa) ! for a 1D column (sounding) of pressures (Pa) and temperatures (K). ! ! - The program tropo_test also corrects for non-sense occuring at high ! latitudes > 60N (60S). ! Theese problems occur at high latitudes on the winter hemisphere because ! of the cold stratosphere and the consequently weak troposphere- ! stratosphere transition of the thermal lapse rate -dT/dz which is ! used to determine the thermal tropopause. ! The problem is much less severe on the winter hemisphere. ! ! Here a fix is implemented applying a factor which scales the ! tropopause pressures at latitudes LAT > 60N (60S) to the ! mean tropopause pressure at 60N (60S), whenever the mean ! tropopause pressure at LAT is lower than that at 60N (60S). ! This idea was formulated in the ECHAM4/CHEM routine xttphwmo by ! D. Nodorp, Ch. Land, B. Steil and R. Hein. ! ! Programmed by Dominik Brunner, ETHZ, Switzerland V.01, 10 Aug 2000 ! Modified by Didier Hauglustaine, IPSL, for LMDZ/INCA, Oct 2000. ! !***************************************************************************** USE AIRPLANE_SRC, ONLY : itrop, ptrop USE INCA_DIM USE MOD_INCA_PARA, ONLY : ij_begin,ij_end, & nbp_para_begin,nbp_para_end,mpi_rank, & plon_omp_begin, plon_omp_end USE MOD_GRID_INCA USE MOD_INCA_OMP_DATA IMPLICIT NONE !---------------------------------------------------------------------- ! ... Dummy args !---------------------------------------------------------------------- REAL, INTENT(IN) :: pmid(PLON_GLO,PLEV) REAL, INTENT(IN) :: temp(PLON_GLO,PLEV) INTEGER, INTENT(IN) :: nx, ny, nz !---------------------------------------------------------------------- ! ... Local variables !---------------------------------------------------------------------- REAL,DIMENSION(nx,ny,nz) :: p3 ! pressures at model layers (full levels) (Pa) REAL,DIMENSION(nx,ny,nz) :: t3 ! temp. at model layer (full levels) (K) INTEGER :: i,j,k,l ! lon, lat and vertical indices REAL, DIMENSION(nx,ny) :: ptropd ! 2D field of tropopause pressures INTEGER, DIMENSION(nx,ny) :: itropd ! 2D field of tropopause level indices REAL, DIMENSION(nx,ny) :: itropdr ! 2D field of tropopause level indices REAL, DIMENSION(PLON_GLO) :: itropr REAL, DIMENSION(nz) :: p1,t1 ! 1D columns of pressures & temperatures INTEGER :: jlat60S, jlat60N ! latitude indices for 60N and 60S REAL, PARAMETER :: rinval=-999.00 ! invalid real number INTEGER :: l300hPa ! index of layer containing 300 hPa level INTEGER :: l120hPa ! index of layer containing 120 hPa level REAL :: ptpm60,ptpmlat,ptpsum ! mean and integrated TP pressures REAL :: rscale ! scaling factor to correct non-sense REAL :: deltalat REAL :: ptrop_glo(PLON_GLO) INTEGER :: itrop_glo(PLON_GLO) INTEGER :: nbbeg_loc,nbend_loc deltalat = 180./FLOAT(ny-1) !T and p to 2D grid CALL grid1dto2d_glo(pmid,p3) CALL grid1dTo2d_glo(temp,t3) !indices of lat=60N and 60S jlat60N=FLOOR(30./deltalat)+1 jlat60S=ny-jlat60N+1 !Calc index of model layer (typically) containing the 300 hPa level !loop from surface to top DO l=1,nz-1 IF (0.5*(LOG(p3(nx/2,ny/2,l))+LOG(p3(nx/2,ny/2,l+1))).LT.LOG(30000.)) & GOTO 100 ENDDO 100 CONTINUE l300hPa=l !Calc index of model layer (typically) containing the 120 hPa level !loop from surface to top DO l=1,nz-1 IF (0.5*(LOG(p3(nx/2,ny/2,l))+LOG(p3(nx/2,ny/2,l+1))).LT.LOG(12000.)) & GOTO 101 ENDDO 101 CONTINUE l120hPa=l !loop over southern hemisphere and get tropopause (loop starts at equator!!) DO j=ny/2,ny ptpsum=0. ! initialize sum of tropopause pressures DO i=1,nx !reverse order of levels because stattrop !requires first index at top of atmosphere DO l=1,nz p1(l)=p3(i,j,nz+1-l) t1(l)=t3(i,j,nz+1-l) ENDDO CALL stattrop(p1,t1,nz,ptropd(i,j),itropd(i,j)) itropd(i,j)=nz-itropd(i,j)+1 ! reverse tropopause level index back !if no tropopause was found (happens only rarely) IF (ptropd(i,j).EQ.rinval) THEN IF (j.GT.jlat60S) THEN ptropd(i,j)=30000. ! set to 300 hPa at high latitudes itropd(i,j)=l300hPa ! index of layer containing 300 hPa level ELSE ptropd(i,j)=12000. ! set to 120 hPa at low latitudes itropd(i,j)=l120hPa ! index of layer containing 120 hPa level END IF ENDIF !sum up tropopause levels to calculate mean ptpsum=ptpsum+ptropd(i,j) ENDDO !Correct nonsense at high latitudes (correction is only !necessary during southern hemisphere winter from Apr until Sep) IF (j.EQ.jlat60S) ptpm60=ptpsum/float(nx) IF (j.GT.jlat60S) THEN ptpmlat=ptpsum/float(nx) !correct ptropd if mean tropopause pressure at this latitude !is lower than at 60S IF (ptpmlat.LT.ptpm60) THEN rscale=ptpm60/ptpmlat DO i=1,nx ptropd(i,j)=ptropd(i,j)*rscale !find corresponding level indices (start loop at surface) DO l=1,nz-1 IF (0.5*(LOG(p3(i,j,l))+LOG(p3(i,j,l+1))) & .LT.LOG(ptropd(i,j))) GOTO 102 ENDDO 102 itropd(i,j)=l ENDDO ENDIF ENDIF ENDDO !loop over northern hemisphere and get tropopause (loop starts at equator!!) DO j=ny/2,1,-1 ptpsum=0. ! initialize sum of tropopause pressures DO i=1,nx !reverse order of levels because stattrop !requires first index at top of atmosphere DO l=1,nz p1(l)=p3(i,j,nz+1-l) t1(l)=t3(i,j,nz+1-l) ENDDO CALL stattrop(p1,t1,nz,ptropd(i,j),itropd(i,j)) itropd(i,j)=nz-itropd(i,j)+1 ! reverse tropopause level index back !if no tropopause was found (happens only rarely) IF (ptropd(i,j).EQ.rinval) THEN IF (j.LT.jlat60N) THEN ptropd(i,j)=30000. ! set to 300 hPa for high latitudes itropd(i,j)=l300hPa ! index of layer containing 300 hPa level ELSE ptropd(i,j)=12000. ! set to 120 hPa for low latitudes itropd(i,j)=l120hPa ! index of layer containing 120 hPa level END IF ENDIF !sum up tropopause levels to calculate mean ptpsum=ptpsum+ptropd(i,j) ENDDO !Correct nonsense at high latitudes (correction is only !necessary during nouthern hemisphere winter from Oct-Mar) IF (j.EQ.jlat60N) ptpm60=ptpsum/float(nx) IF (j.LT.jlat60N) THEN ptpmlat=ptpsum/float(nx) !correct ptropd if mean tropopause pressure at this latitude !is lower than at 60S IF (ptpmlat.LT.ptpm60) THEN rscale=ptpm60/ptpmlat DO i=1,nx ptropd(i,j)=ptropd(i,j)*rscale !find corresponding level indices (start loop at surface) DO l=1,nz-1 IF (0.5*(LOG(p3(i,j,l))+LOG(p3(i,j,l+1))) & .LT.LOG(ptropd(i,j))) GOTO 103 ENDDO 103 itropd(i,j)=l ENDDO ENDIF ENDIF ENDDO ! borne pour le passage GLO -> LOC nbbeg_loc = nbp_para_begin(mpi_rank)+plon_omp_begin-1 nbend_loc = nbp_para_begin(mpi_rank)+plon_omp_end-1 !Back to physiq grid 1D CALL grid2dto1d_glo(ptropd,ptrop_glo) ptrop=ptrop_glo(nbbeg_loc:nbend_loc) itropdr = FLOAT(itropd) CALL grid2dto1d_glo(itropdr,itropr) itrop_glo = INT(itropr) itrop=itrop_glo(nbbeg_loc:nbend_loc) END SUBROUTINE FDTROPOPAUSE SUBROUTINE PINTERP(ncol, & fieldin, & pin, & nin, & fieldout, & pout, & nout, & index, & entered) !----------------------------------------------------------------------- ! ... Interpolates on model pressure levels ! Stacy Walters and Didier Hauglustaine, 1999. !----------------------------------------------------------------------- IMPLICIT NONE !----------------------------------------------------------------------- ! ... Dummy args !----------------------------------------------------------------------- LOGICAL, INTENT(in) :: entered INTEGER, INTENT(in) :: nin INTEGER, INTENT(in) :: nout INTEGER, INTENT(in) :: ncol INTEGER, INTENT(inout) :: INDEX(ncol,nout) REAL, INTENT(in) :: pin(ncol,nin) REAL, INTENT(in) :: fieldin(ncol,nin) REAL, INTENT(in) :: pout(ncol,nout) REAL, INTENT(out) :: fieldout(ncol,nout) !----------------------------------------------------------------------- ! ... Local variables !----------------------------------------------------------------------- INTEGER :: i, j, k REAL :: deltp(ncol,nout) REAL :: pinratio(ncol,nin) IF ( .NOT. entered) THEN DO i = 1, ncol DO k = 1, nout DO j = nin-1, 1, -1 IF (pout(i,k) <= pin(i,j)) THEN IF (pout(i,k) > pin(i,j+1)) THEN index(i,k) = j ENDIF ENDIF ENDDO ENDDO ENDDO ENDIF DO i = 1, ncol DO j = 1, nin-1 pinratio(i,j) = 1./LOG(pin(i,j)/pin(i,j+1)) END DO END DO DO i = 1, ncol DO k = 1, nout IF (pout(i,k) <= pin(i,nin)) THEN fieldout(i,k) = fieldin(i,nin) ELSE IF (pout(i,k) >= pin(i,1)) THEN fieldout(i,k) = fieldin(i,1) ELSE deltp(i,k) = LOG(pout(i,k)/pin(i,index(i,k)+1)) * pinratio(i,index(i,k)) fieldout(i,k) = fieldin(i,index(i,k)+1) & + deltp(i,k) * (fieldin(i,index(i,k))-fieldin(i,index(i,k)+1)) ENDIF ENDDO ENDDO END SUBROUTINE PINTERP SUBROUTINE stattrop(pfull, tfull, nlev, ptropd, itropd) !********************************************************************************** ! ! input : p, t real(nlev) ! input : nlev real ! output: ptropd real ! output: itropd integer ! ! programmed by Dominik Brunner V1.0 Aug 2000 ! ! built upon routine stattrop by Peter van Velthoven, ! KNMI, The Netherlands ! and on the ECHAM4/CHEM routine xttphwmo by ! Thomas Reichle, Christine Land, B. Steil and R. Hein, DLR ! ! purpose ! ------- ! stattrop computes the pressure (Pa) at the thermal (static) tropopause ! (TP) for a 1D column (sounding) of pressures and temperatures following ! the definition of the height of the tropopause as postulated by WMO (1957). ! ! ATTENTION: In the current formulation of the program the first ! level (index 1) must be at the top of the atmosphere and the ! last level (index nlev) at the surface ! ! interface ! --------- ! call stattrop(pfull, tfull, nlev, ptropd, itropd) ! - Input ! nlev : number of vertical levels ! pfull: pressure in 1D column at nlev full levels (layers) ! tfull: temperature in 1D column at nlev full levels ! - Output ! ptropd: height of the tropopause in Pa ! itropd: index of layer containing the tropopause ! ! methods ! ------- ! - Lapse rate gamma = -dT/dz ! Using the hydrostatic approximation ! ! dz = -dp/(g*rho) = -dp/p * R/g * T = -dlnp * R/g * T ! ! we get -dT/dz = dT/T * g/R *1/dlnp = dlnT/dlnp ! ! - Variables are assumed to vary linearly in log(pressure) ! - The tropopause is the lowest level above 450 hPa fullfilling ! the WMO criterium. If ptropd is < 85 hPa it is set to 85 hPa. ! If no tropopause is found ptropd is set to -999. ! ! references ! ---------- ! ! - WMO (1992): International meteorological vocabulary, Genf, 784pp.: ! ! 1. The first tropopause is defined as the lowest level at which ! the lapse rate decreases to 2 deg C per kilometer or less, ! provided also the average lapse rate between this level and ! all higher levels within 2 kilometers does not exceed 2 deg C ! ! - Randel WJ, Wu F, Gaffen DJ, Interannual variability of the tropical ! tropopause derived from radiosonde data and NCEP reanalyses, ! JOURNAL OF GEOPHYSICAL RESEARCH, 105, 15509-15523, 2000. ! ! The following webpage clearifies the calculation of the tropopause ! in the NCEP reanalysis: http://dss.ucar.edu/pub/reanalysis/FAQ.html ! ! For comparison NCEP reanalysis tropopause pressures can be obtained ! from http://www.cdc.noaa.gov/cdc/reanalysis/reanalysis.shtml ! !********************************************************************************** USE INCA_DIM IMPLICIT NONE INTEGER nlev REAL pfull(nlev), tfull(nlev) REAL ptropd, p2km, lapseavg, lapsesum INTEGER ilev, ilev1, itropd, count, l REAL lapse(PLEV), phalf(PLEV) REAL grav, rgas, const, lapsec, rinval PARAMETER (grav=9.80665, rgas=287.04, & const=1000.*grav/rgas, & lapsec=2.0, rinval=-999.0 ) ! min. and max. allowed tropopause pressure REAL pmin, pmax PARAMETER (pmin=8500.,pmax=45000.) ! deltaz is layer thickness used for second tropopause criterium REAL deltaz PARAMETER (deltaz=2.) LOGICAL found REAL p1, p2 , weight ptropd = rinval itropd = 0 found = .FALSE. ! Loop from model top to surface and calculate lapse rate gamma=-dT/dz (K/km) DO ilev = 1, nlev-1 ilev1 = ilev + 1 lapse(ilev) = LOG(tfull(ilev))- LOG(tfull(ilev1)) lapse(ilev) = lapse(ilev) / & (LOG(pfull(ilev))- LOG(pfull(ilev1))) lapse(ilev) = const * lapse(ilev) phalf(ilev) = (pfull(ilev) + pfull(ilev1))*0.5 END DO ! Loop from surface to top to find (lowest) tropopause DO ilev = nlev-2, 1, -1 ! **** 1st test: -dT/dz less than 2 K/km and pressure LT pmax? **** !Modified 07/2001 -- D. Hauglustaine IF (lapse(ilev).LT.lapsec.AND.pfull(ilev).LT.pmax.AND.ilev.GE.5) THEN IF (.NOT.found) THEN ! Interpolate tropopause pressure linearly in log(p) p1 = LOG(phalf(ilev)) p2 = LOG(phalf(ilev+1)) !Modified 09/2001 -- L. Jourdain IF ( (lapse(ilev).NE.lapse(ilev+1)) .AND. (lapse(ilev+1).GE.lapsec) ) THEN weight = (lapsec - lapse(ilev)) / & (lapse(ilev+1)-lapse(ilev)) ptropd = EXP(p1 + weight * (p2-p1)) ! tropopause pressure ELSE ptropd = phalf(ilev) END IF ! *** 2nd test: average -dT/dz in layer deltaz above TP must not be greater than 2K/km p2km = EXP(LOG(ptropd)-deltaz*const/tfull(ilev)) ! press 2 km above TP lapsesum = 0.0 ! init avg. lapse rate 2 km above TP lapseavg=1.e9 count = 0 DO l=ilev,1,-1 IF (phalf(l).LT.p2km) GOTO 100 lapsesum=lapsesum+lapse(l) count=count+1 lapseavg=lapsesum/float(count) END DO 100 CONTINUE found=lapseavg.LT.lapsec ! Discard tropopause? IF (.NOT.found) THEN ptropd=rinval ELSE ptropd=MAX(ptropd,pmin) ! ptropd must be GE 85 hPa itropd=ilev ! assign level index END IF END IF END IF END DO RETURN END SUBROUTINE stattrop subroutine CALC_PV(latitude_deg,paprs,pplay,ts,t,rot) ! This subroutine consists in deriving Ertel's potential vorticity (PV) in ! potential vorticity units (PVU), i.e. 1 PVU = 1E-6 K.m^2/kg/s. ! It accounts for the vertical potential temperature gradient (=static stability) ! and the (almost) horizontal wind vorticity. ! The equation and clear explanations are available in Gettelman et al. (2011, Rev. Geophys), ! a review paper on the extratropical UTLS. ! Added by Yann Cohen, November 2019. ! Created by Yann Cohen USE CONST_MOD USE XIOS_INCA IMPLICIT NONE !------------------------------------------------------------------------------- ! Arguments: REAL, INTENT(IN) :: latitude_deg(PLON) ! latitude REAL, INTENT(IN) :: paprs(PLON,PLEVP) !--- Cells-edges pressure REAL, INTENT(IN) :: pplay(PLON,PLEV) !--- Cells-centers pressure REAL, INTENT(IN) :: ts(PLON) !--- Surface temperature REAL, INTENT(IN) :: t(PLON,PLEV) !--- Cells-centers temperature REAL, INTENT(IN) :: rot(PLON,PLEV) !--- Cells-centers relative vorticity !------------------------------------------------------------------------------- ! Local variables: include "YOMCST_I.h" INTEGER :: i, k REAL :: al,al2top REAL,PARAMETER :: preff=101325 ! Unit: Pa. REAL, DIMENSION(PLON,PLEVP):: tpot_edg REAL, DIMENSION(PLON,PLEV) :: tpot_cen REAL, DIMENSION(PLON,PLEV) :: PV !------------------------------------------------------------------------------- PV(:,:) = 0. !--- POTENTIAL TEMPERATURE AT CELLS CENTERS AND INTERFACES DO i = 1,PLON tpot_cen(i,1) = t(i,1)*(preff/pplay(i,1))**RKAPPA tpot_edg(i,1) = ts(i) *(preff/paprs(i,1))**RKAPPA DO k=2,PLEV tpot_cen(i,k) = t(i,k)*(preff/pplay(i,k))**RKAPPA al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k)) tpot_edg(i,k) = (t(i,k-1)+al*(t(i,k)-t(i,k-1)))*(preff/paprs(i,k))**RKAPPA END DO ! Last ingredient for the PV field on the whole vertical grid: ! linear extrapolation up to the summit, for cell-edges potential temperature. ! The subsequent Theta vertical gradient is then used for the PV calculation at the top (full) level. ! al2top=LOG(paprs(i,PLEV)/paprs(i,PLEVP))/LOG(paprs(i,PLEV)/pplay(i,PLEV)) ! tpot_edg(i,PLEVP) = tpot_edg(i,PLEV) + al2top*(tpot_cen(i,PLEV)-tpot_edg(i,PLEV)) END DO !--- ERTEL POTENTIAL VORTICITY AT CELLS CENTERS (except in top layer) DO i = 1, PLON DO k= 1, PLEV-1 PV(i,k)=-1.E6*RG*(rot(i,k)+2.*ROMEGA*SIN(latitude_deg(i)*RPI/180.))& * (tpot_edg(i,k+1)-tpot_edg(i,k)) / (paprs(i,k+1)-paprs(i,k)) END DO END DO CALL xios_inca_change_context("inca") CALL xios_inca_send_field("PV", PV) CALL xios_inca_change_context("LMDZ") end subroutine CALC_PV