[6610] | 1 | !$Id: set_ub_vals.F90 157 2010-01-18 14:21:59Z acosce $ |
---|
| 2 | !! ========================================================================= |
---|
| 3 | !! INCA - INteraction with Chemistry and Aerosols |
---|
| 4 | !! |
---|
| 5 | !! Copyright Laboratoire des Sciences du Climat et de l'Environnement (LSCE) |
---|
| 6 | !! Unite mixte CEA-CNRS-UVSQ |
---|
| 7 | !! |
---|
| 8 | !! Contributors to this INCA subroutine: |
---|
| 9 | !! |
---|
| 10 | !! Didier Hauglustaine, LSCE, hauglustaine@cea.fr |
---|
| 11 | !! |
---|
| 12 | !! Anne Cozic, LSCE, anne.cozic@cea.fr |
---|
| 13 | !! Yann Meurdesoif, LSCE, yann.meurdesoif@cea.fr |
---|
| 14 | !! |
---|
| 15 | !! This software is a computer program whose purpose is to simulate the |
---|
| 16 | !! atmospheric gas phase and aerosol composition. The model is designed to be |
---|
| 17 | !! used within a transport model or a general circulation model. This version |
---|
| 18 | !! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts |
---|
| 19 | !! for emissions, transport (resolved and sub-grid scale), photochemical |
---|
| 20 | !! transformations, and scavenging (dry deposition and washout) of chemical |
---|
| 21 | !! species and aerosols interactively in the GCM. Several versions of the INCA |
---|
| 22 | !! model are currently used depending on the envisaged applications with the |
---|
| 23 | !! chemistry-climate model. |
---|
| 24 | !! |
---|
| 25 | !! This software is governed by the CeCILL license under French law and |
---|
| 26 | !! abiding by the rules of distribution of free software. You can use, |
---|
| 27 | !! modify and/ or redistribute the software under the terms of the CeCILL |
---|
| 28 | !! license as circulated by CEA, CNRS and INRIA at the following URL |
---|
| 29 | !! "http://www.cecill.info". |
---|
| 30 | !! |
---|
| 31 | !! As a counterpart to the access to the source code and rights to copy, |
---|
| 32 | !! modify and redistribute granted by the license, users are provided only |
---|
| 33 | !! with a limited warranty and the software's author, the holder of the |
---|
| 34 | !! economic rights, and the successive licensors have only limited |
---|
| 35 | !! liability. |
---|
| 36 | !! |
---|
| 37 | !! In this respect, the user's attention is drawn to the risks associated |
---|
| 38 | !! with loading, using, modifying and/or developing or reproducing the |
---|
| 39 | !! software by the user in light of its specific status of free software, |
---|
| 40 | !! that may mean that it is complicated to manipulate, and that also |
---|
| 41 | !! therefore means that it is reserved for developers and experienced |
---|
| 42 | !! professionals having in-depth computer knowledge. Users are therefore |
---|
| 43 | !! encouraged to load and test the software's suitability as regards their |
---|
| 44 | !! requirements in conditions enabling the security of their systems and/or |
---|
| 45 | !! data to be ensured and, more generally, to use and operate it in the |
---|
| 46 | !! same conditions as regards security. |
---|
| 47 | !! |
---|
| 48 | !! The fact that you are presently reading this means that you have had |
---|
| 49 | !! knowledge of the CeCILL license and that you accept its terms. |
---|
| 50 | !! ========================================================================= |
---|
| 51 | |
---|
| 52 | #include <inca_define.h> |
---|
| 53 | |
---|
| 54 | #if defined(AERONLY) || defined(GES) |
---|
| 55 | |
---|
| 56 | |
---|
| 57 | |
---|
| 58 | |
---|
| 59 | SUBROUTINE XIOS_OXYDANT_READ(calday) |
---|
| 60 | !---------------------------------------------------------------------- |
---|
| 61 | ! ... Read OXYDANT distributions |
---|
| 62 | ! Didier Hauglustaine, IPSL, 1999. |
---|
| 63 | !---------------------------------------------------------------------- |
---|
| 64 | |
---|
| 65 | USE OXYDANT_COM |
---|
| 66 | USE MOD_GRID_INCA |
---|
| 67 | USE INCA_DIM |
---|
| 68 | USE MOD_INCA_PARA |
---|
| 69 | USE PRINT_INCA |
---|
| 70 | USE XIOS_INCA |
---|
| 71 | IMPLICIT NONE |
---|
| 72 | |
---|
| 73 | INCLUDE 'netcdf.inc' |
---|
| 74 | |
---|
| 75 | !---------------------------------------------------------------------- |
---|
| 76 | ! ... Dummy args |
---|
| 77 | !---------------------------------------------------------------------- |
---|
| 78 | REAL, INTENT(IN) :: calday |
---|
| 79 | |
---|
| 80 | !---------------------------------------------------------------------- |
---|
| 81 | ! ... Local variables |
---|
| 82 | !---------------------------------------------------------------------- |
---|
| 83 | INTEGER :: iret |
---|
| 84 | INTEGER :: varid |
---|
| 85 | INTEGER :: oldlotime, oldhitime |
---|
| 86 | INTEGER, DIMENSION(4) :: start3, count3 |
---|
| 87 | LOGICAL, SAVE :: first = .TRUE. |
---|
| 88 | !$OMP THREADPRIVATE(first) |
---|
| 89 | |
---|
| 90 | REAL, DIMENSION(12) :: timeo_inter |
---|
| 91 | REAL :: dt, dt1, tint |
---|
| 92 | INTEGER :: lonlev |
---|
| 93 | |
---|
| 94 | |
---|
| 95 | IF (first) THEN |
---|
| 96 | |
---|
| 97 | |
---|
| 98 | |
---|
| 99 | ALLOCATE(timeo(ntime_oxyd)) |
---|
| 100 | ALLOCATE(o3oxyd_year(PLON, PLEV, ntime_oxyd)) |
---|
| 101 | ALLOCATE(o1doxyd_year(PLON, PLEV, ntime_oxyd)) |
---|
| 102 | ALLOCATE(hno3oxyd_year(PLON, PLEV, ntime_oxyd)) |
---|
| 103 | ALLOCATE(ohoxyd_year(PLON, PLEV, ntime_oxyd)) |
---|
| 104 | ALLOCATE(h2o2oxyd_year(PLON, PLEV, ntime_oxyd)) |
---|
| 105 | ALLOCATE(no2oxyd_year(PLON, PLEV, ntime_oxyd)) |
---|
| 106 | ALLOCATE(no3oxyd_year(PLON, PLEV, ntime_oxyd)) |
---|
| 107 | |
---|
| 108 | ! write(lunout,*) '--------------' |
---|
| 109 | ! write(lunout,*) 'lecture oxyd : ' |
---|
| 110 | ! write(lunout,*) 'PLON = ', PLON |
---|
| 111 | ! write(lunout,*) 'PLEV = ', PLEV |
---|
| 112 | ! write(lunout,*) 'ntime_oxy = ', ntime_oxyd |
---|
| 113 | ! write(lunout,*) '--------------' |
---|
| 114 | |
---|
| 115 | CALL xios_inca_change_context("inca") |
---|
| 116 | |
---|
| 117 | CALL xios_inca_recv_field_glo("oxydtime_read", timeo ) |
---|
| 118 | |
---|
| 119 | CALL xios_inca_recv_field("O3_interp" , o3oxyd_year ) |
---|
| 120 | CALL xios_inca_recv_field("O1D_interp" , o1doxyd_year ) |
---|
| 121 | CALL xios_inca_recv_field("HNO3_interp" , hno3oxyd_year) |
---|
| 122 | CALL xios_inca_recv_field("OH_interp" , ohoxyd_year ) |
---|
| 123 | CALL xios_inca_recv_field("H2O2_interp" , h2o2oxyd_year) |
---|
| 124 | CALL xios_inca_recv_field("NO2_interp" , no2oxyd_year ) |
---|
| 125 | CALL xios_inca_recv_field("NO3_interp" , no3oxyd_year ) |
---|
| 126 | |
---|
| 127 | CALL xios_inca_send_field("NO3_oxyd_read", no3oxyd_year(:,:,12)) |
---|
| 128 | |
---|
| 129 | CALL xios_inca_change_context("LMDZ") |
---|
| 130 | |
---|
| 131 | |
---|
| 132 | |
---|
| 133 | first = .false. |
---|
| 134 | ENDIF |
---|
| 135 | |
---|
| 136 | |
---|
| 137 | lonlev = PLON * nlevo |
---|
| 138 | |
---|
| 139 | |
---|
| 140 | ! ... Check to see if model time still bounded by data set |
---|
| 141 | oldlotime = lotime |
---|
| 142 | oldhitime = hitime |
---|
| 143 | timeo_inter(:) = mod(timeo(:)/86400, 365.) |
---|
| 144 | |
---|
| 145 | |
---|
| 146 | CALL FINDPLB (timeo_inter, ntime_oxyd, calday, lotime) |
---|
| 147 | |
---|
| 148 | IF ( cyclical ) THEN |
---|
| 149 | hitime = MOD(lotime,ntime_oxyd) + 1 |
---|
| 150 | ELSE |
---|
| 151 | hitime = lotime + 1 |
---|
| 152 | END IF |
---|
| 153 | ! -------------------- interpolation sur le temps |
---|
| 154 | !---------------------------------------------------------------------- |
---|
| 155 | ! Note: 365 day year inconsistent with LMDz (360 days) !!! |
---|
| 156 | !---------------------------------------------------------------------- |
---|
| 157 | ! ... First interpolate linearly on current model time |
---|
| 158 | ! |
---|
| 159 | ! |
---|
| 160 | IF ( timeo(hitime) < timeo(lotime) ) THEN |
---|
| 161 | dt = 365. - mod(timeo(lotime)/(86400),365.) + mod(timeo(hitime)/(86400),365.) |
---|
| 162 | IF (calday <= mod(timeo(hitime)/(86400),365.) ) THEN |
---|
| 163 | dt1 = 365. - mod(timeo(lotime)/(86400),365.) + calday |
---|
| 164 | ELSE |
---|
| 165 | dt1 = calday - mod(timeo(lotime)/(86400),365.) |
---|
| 166 | END IF |
---|
| 167 | ELSE |
---|
| 168 | dt = mod(timeo(hitime)/(86400),365.) - mod(timeo(lotime)/(86400),365.) |
---|
| 169 | dt1 = calday - mod(timeo(lotime)/(86400),365.) |
---|
| 170 | END IF |
---|
| 171 | |
---|
| 172 | tint = dt1/dt |
---|
| 173 | |
---|
| 174 | ohoxyd(:,:) = ohoxyd_year(:,:,lotime) + (ohoxyd_year(:,:,hitime)-ohoxyd_year(:,:,lotime))*tint |
---|
| 175 | o1doxyd(:,:) = o1doxyd_year(:,:,lotime) + (o1doxyd_year(:,:,hitime)-o1doxyd_year(:,:,lotime))*tint |
---|
| 176 | o3oxyd(:,:) = o3oxyd_year(:,:,lotime) + (o3oxyd_year(:,:,hitime)-o3oxyd_year(:,:,lotime))*tint |
---|
| 177 | no3oxyd(:,:) = no3oxyd_year(:,:,lotime) + (no3oxyd_year(:,:,hitime)-no3oxyd_year(:,:,lotime))*tint |
---|
| 178 | h2o2oxyd(:,:) = h2o2oxyd_year(:,:,lotime) + (h2o2oxyd_year(:,:,hitime)-h2o2oxyd_year(:,:,lotime))*tint |
---|
| 179 | hno3oxyd(:,:) = hno3oxyd_year(:,:,lotime) + (hno3oxyd_year(:,:,hitime)-hno3oxyd_year(:,:,lotime))*tint |
---|
| 180 | no2oxyd(:,:) = no2oxyd_year(:,:,lotime) + (no2oxyd_year(:,:,hitime)-no2oxyd_year(:,:,lotime))*tint |
---|
| 181 | |
---|
| 182 | CALL xios_inca_change_context("inca") |
---|
| 183 | CALL xios_inca_send_field("O3_oxyd" , o3oxyd ) |
---|
| 184 | CALL xios_inca_send_field("O1D_oxyd" , o1doxyd ) |
---|
| 185 | CALL xios_inca_send_field("HNO3_oxyd" , hno3oxyd ) |
---|
| 186 | CALL xios_inca_send_field("OH_oxyd" , ohoxyd ) |
---|
| 187 | CALL xios_inca_send_field("H2O2_oxyd" , h2o2oxyd ) |
---|
| 188 | CALL xios_inca_send_field("NO2_oxyd" , no2oxyd ) |
---|
| 189 | CALL xios_inca_send_field("NO3_oxyd" , no3oxyd ) |
---|
| 190 | CALL xios_inca_change_context("LMDZ") |
---|
| 191 | |
---|
| 192 | |
---|
| 193 | END SUBROUTINE XIOS_OXYDANT_READ |
---|
| 194 | |
---|
| 195 | |
---|
| 196 | #else |
---|
| 197 | ! IF not defined aeronly and ges |
---|
| 198 | |
---|
| 199 | SUBROUTINE OZCLIM_READ(calday) |
---|
| 200 | !---------------------------------------------------------------------- |
---|
| 201 | ! ... Read ozone distributions |
---|
| 202 | ! Didier Hauglustaine, IPSL, 1999. |
---|
| 203 | !---------------------------------------------------------------------- |
---|
| 204 | |
---|
| 205 | USE O3CLIM_COM |
---|
| 206 | USE INCA_DIM |
---|
| 207 | USE MOD_INCA_PARA |
---|
| 208 | USE PRINT_INCA |
---|
| 209 | |
---|
| 210 | IMPLICIT NONE |
---|
| 211 | |
---|
| 212 | INCLUDE 'netcdf.inc' |
---|
| 213 | |
---|
| 214 | !---------------------------------------------------------------------- |
---|
| 215 | ! ... Dummy args |
---|
| 216 | !---------------------------------------------------------------------- |
---|
| 217 | REAL, INTENT(IN) :: calday |
---|
| 218 | |
---|
| 219 | !---------------------------------------------------------------------- |
---|
| 220 | ! ... Local variables |
---|
| 221 | !---------------------------------------------------------------------- |
---|
| 222 | INTEGER :: iret |
---|
| 223 | INTEGER :: varid |
---|
| 224 | INTEGER :: oldlotime, oldhitime |
---|
| 225 | INTEGER, DIMENSION(3) :: start, count |
---|
| 226 | REAL :: o3climbd_glo(nbp_glo,nlevc,2) |
---|
| 227 | |
---|
| 228 | ! ... Check to see if model time still bounded by data set |
---|
| 229 | oldlotime = lotime |
---|
| 230 | oldhitime = hitime |
---|
| 231 | CALL FINDPLB (timec, ntimec, calday, lotime) |
---|
| 232 | IF ( cyclical ) THEN |
---|
| 233 | hitime = MOD(lotime,ntimec) + 1 |
---|
| 234 | ELSE |
---|
| 235 | hitime = lotime + 1 |
---|
| 236 | END IF |
---|
| 237 | |
---|
| 238 | ! ... Read new data if needed |
---|
| 239 | IF ( hitime /= oldhitime ) THEN |
---|
| 240 | loind = hiind |
---|
| 241 | hiind = MOD( loind, 2) + 1 |
---|
| 242 | |
---|
| 243 | !$OMP MASTER |
---|
| 244 | IF (is_mpi_root) THEN |
---|
| 245 | iret = nf_inq_varid(ncidc, 'O3', varid) |
---|
| 246 | CALL check_err(iret, 'OZCLIM_READ','problem to check varid O3') |
---|
| 247 | |
---|
| 248 | count(1) = nbp_glo |
---|
| 249 | count(2) = nlevc |
---|
| 250 | count(3) = 1 |
---|
| 251 | start(1) = 1 |
---|
| 252 | start(2) = 1 |
---|
| 253 | |
---|
| 254 | WRITE(lunout,*) 'OZCLIM_READ : read new data for time ', timec(hitime) |
---|
| 255 | start(3) = hitime |
---|
| 256 | iret = nf_get_vara_double(ncidc, varid, start, count,o3climbd_glo(1,1,hiind)) |
---|
| 257 | CALL check_err(iret, 'OZCLIM_READ', 'problem to read O3 hiind') |
---|
| 258 | ENDIF |
---|
| 259 | !$OMP END MASTER |
---|
| 260 | CALL scatter(o3climbd_glo(:,:,hiind),o3climbd(:,:,hiind)) |
---|
| 261 | |
---|
| 262 | |
---|
| 263 | IF ( lotime /= oldhitime ) THEN |
---|
| 264 | !$OMP MASTER |
---|
| 265 | IF (is_mpi_root) THEN |
---|
| 266 | WRITE(lunout,*) 'OZCLIM_READ : read new data for time ', timec(lotime) |
---|
| 267 | start(3) = lotime |
---|
| 268 | iret = nf_get_vara_double(ncidc, varid, start, count, o3climbd_glo(1,1,loind)) |
---|
| 269 | CALL check_err(iret, 'OZCLIM_READ', 'problem to read O3 loind') |
---|
| 270 | |
---|
| 271 | ENDIF |
---|
| 272 | !$OMP END MASTER |
---|
| 273 | CALL scatter(o3climbd_glo(:,:,loind),o3climbd(:,:,loind)) |
---|
| 274 | |
---|
| 275 | END IF |
---|
| 276 | |
---|
| 277 | END IF |
---|
| 278 | |
---|
| 279 | END SUBROUTINE OZCLIM_READ |
---|
| 280 | |
---|
| 281 | SUBROUTINE OZCLIM_INTERP(calday, pmid) |
---|
| 282 | !---------------------------------------------------------------------- |
---|
| 283 | ! ... Interpolate ozone on model time and on pressure levels |
---|
| 284 | ! Didier Hauglustaine, IPSL, 1999. |
---|
| 285 | !---------------------------------------------------------------------- |
---|
| 286 | |
---|
| 287 | USE O3CLIM_COM |
---|
| 288 | USE INCA_DIM |
---|
| 289 | IMPLICIT NONE |
---|
| 290 | |
---|
| 291 | !---------------------------------------------------------------------- |
---|
| 292 | ! ... Dummy args |
---|
| 293 | !---------------------------------------------------------------------- |
---|
| 294 | REAL, INTENT(IN) :: calday |
---|
| 295 | REAL, INTENT(IN) :: pmid(PLON,PLEV) |
---|
| 296 | |
---|
| 297 | !---------------------------------------------------------------------- |
---|
| 298 | ! ... Local variables |
---|
| 299 | !---------------------------------------------------------------------- |
---|
| 300 | REAL :: dt, dt1, tint |
---|
| 301 | REAL :: rma = 48./28. |
---|
| 302 | REAL :: ploc(klonc,nlevc) |
---|
| 303 | INTEGER :: lonlev |
---|
| 304 | INTEGER :: i |
---|
| 305 | INTEGER :: index(PLON,PLEV) |
---|
| 306 | |
---|
| 307 | index = -9999 |
---|
| 308 | lonlev = klonc * nlevc |
---|
| 309 | |
---|
| 310 | DO i = 1, klonc |
---|
| 311 | ploc(i,:) = presc(:) |
---|
| 312 | END DO |
---|
| 313 | |
---|
| 314 | !---------------------------------------------------------------------- |
---|
| 315 | ! Note: 365 day year inconsistent with LMDz (360 days) !!! |
---|
| 316 | !---------------------------------------------------------------------- |
---|
| 317 | ! ... First interpolate linearly on current model time |
---|
| 318 | |
---|
| 319 | IF ( timec(hitime) < timec(lotime) ) THEN |
---|
| 320 | dt = 365. - timec(lotime) + timec(hitime) |
---|
| 321 | IF (calday <= timec(hitime) ) THEN |
---|
| 322 | dt1 = 365. - timec(lotime) + calday |
---|
| 323 | ELSE |
---|
| 324 | dt1 = calday - timec(lotime) |
---|
| 325 | END IF |
---|
| 326 | ELSE |
---|
| 327 | dt = timec(hitime) - timec(lotime) |
---|
| 328 | dt1 = calday - timec(lotime) |
---|
| 329 | END IF |
---|
| 330 | |
---|
| 331 | tint = dt1/dt |
---|
| 332 | |
---|
| 333 | CALL linintp (lonlev, & |
---|
| 334 | 0., & |
---|
| 335 | 1., & |
---|
| 336 | tint, & |
---|
| 337 | o3climbd(1,1,loind), & |
---|
| 338 | o3climbd(1,1,hiind), & |
---|
| 339 | o3climcr) |
---|
| 340 | |
---|
| 341 | ! ... Put on model pressure levels |
---|
| 342 | |
---|
| 343 | CALL pinterp (PLON, & |
---|
| 344 | o3climcr, & |
---|
| 345 | ploc, & |
---|
| 346 | nlevc, & |
---|
| 347 | o3clim, & |
---|
| 348 | pmid, & |
---|
| 349 | PLEV, & |
---|
| 350 | index, & |
---|
| 351 | .false. ) |
---|
| 352 | |
---|
| 353 | ! ... vmr to mmr |
---|
| 354 | |
---|
| 355 | o3clim = o3clim * rma |
---|
| 356 | |
---|
| 357 | END SUBROUTINE OZCLIM_INTERP |
---|
| 358 | |
---|
| 359 | |
---|
| 360 | SUBROUTINE OZRESET (mmr, pmid, temp) |
---|
| 361 | !---------------------------------------------------------------------- |
---|
| 362 | ! ... Reset O3 boundary conditions for use in radiation |
---|
| 363 | ! Didier Hauglustaine, IPSL, 2002. |
---|
| 364 | !---------------------------------------------------------------------- |
---|
| 365 | |
---|
| 366 | USE O3CLIM_COM |
---|
| 367 | USE SPECIES_NAMES |
---|
| 368 | USE OXYDANT_COM |
---|
| 369 | USE AIRPLANE_SRC, ONLY : itrop |
---|
| 370 | USE CHEM_MODS, ONLY : adv_mass |
---|
| 371 | USE INCA_DIM |
---|
| 372 | IMPLICIT NONE |
---|
| 373 | |
---|
| 374 | !---------------------------------------------------------------------- |
---|
| 375 | ! ... Dummy args |
---|
| 376 | !---------------------------------------------------------------------- |
---|
| 377 | REAL, INTENT(IN) :: pmid(PLON,PLEV) |
---|
| 378 | REAL, INTENT(IN) :: temp(PLON,PLEV) |
---|
| 379 | REAL, INTENT(INOUT) :: mmr(PLON,PLEV,PCNST) |
---|
| 380 | |
---|
| 381 | !---------------------------------------------------------------------- |
---|
| 382 | ! ... Local variables |
---|
| 383 | !---------------------------------------------------------------------- |
---|
| 384 | INTEGER :: i, k |
---|
| 385 | INTEGER :: ktrop, krelax, k380 |
---|
| 386 | REAL, PARAMETER :: dry_mass = 28.966 |
---|
| 387 | REAL :: theta(PLON,PLEV) |
---|
| 388 | |
---|
| 389 | theta(:,:) = temp(:,:)*(1.e5/pmid(:,:))**0.286 |
---|
| 390 | |
---|
| 391 | DO i = 1, PLON |
---|
| 392 | |
---|
| 393 | ktrop = MAX(1,itrop(i)) |
---|
| 394 | |
---|
| 395 | DO k = ktrop, PLEV |
---|
| 396 | IF (theta(i,k) >= 380.) THEN |
---|
| 397 | k380 = k |
---|
| 398 | EXIT |
---|
| 399 | ENDIF |
---|
| 400 | ENDDO |
---|
| 401 | |
---|
| 402 | ! krelax = MIN(ktrop+2,PLEV) ! 2 levels above tropopause |
---|
| 403 | krelax = k380 ! 380K theta surface |
---|
| 404 | |
---|
| 405 | mmr(i,krelax:PLEV,id_O3) = o3clim(i,krelax:PLEV) |
---|
| 406 | END DO |
---|
| 407 | |
---|
| 408 | END SUBROUTINE OZRESET |
---|
| 409 | |
---|
| 410 | SUBROUTINE FIELD_PREP (mmr) |
---|
| 411 | !---------------------------------------------------------------------- |
---|
| 412 | ! ... Prepare fields to be used in LMDz radiation |
---|
| 413 | ! Didier Hauglustaine, IPSL, 2002. 2018. |
---|
| 414 | !---------------------------------------------------------------------- |
---|
| 415 | |
---|
| 416 | USE SPECIES_NAMES |
---|
| 417 | USE CHEM_MODS, ONLY : o3_inca, h2o_inca, ch4_inca, n2o_inca, cfc11_inca, cfc12_inca, & |
---|
| 418 | ch4_inca_2d, n2o_inca_2d, cfc11_inca_2d, cfc12_inca_2d, adv_mass |
---|
| 419 | USE INCA_DIM |
---|
| 420 | IMPLICIT NONE |
---|
| 421 | |
---|
| 422 | !---------------------------------------------------------------------- |
---|
| 423 | ! ... Dummy args |
---|
| 424 | !---------------------------------------------------------------------- |
---|
| 425 | REAL, INTENT(IN) :: mmr(PLON,PLEV,PCNST) |
---|
| 426 | |
---|
| 427 | !---------------------------------------------------------------------- |
---|
| 428 | ! ... Local variables |
---|
| 429 | !---------------------------------------------------------------------- |
---|
| 430 | INTEGER :: jj, jk, jki, jkf |
---|
| 431 | INTEGER, PARAMETER :: ng1 = 2, ng1p1 = ng1 + 1 !From raddimlw.h |
---|
| 432 | |
---|
| 433 | ! Fields to be used in LMDz in mass mixing ratio |
---|
| 434 | #ifdef NMHC |
---|
| 435 | o3_inca(:,:) = mmr(:,:,id_O3) |
---|
| 436 | #endif |
---|
| 437 | ! ch4_inca(:,:) = mmr(:,:,id_CH4) |
---|
| 438 | ! n2o_inca(:,:) = mmr(:,:,id_N2O) |
---|
| 439 | ! cfc11_inca(:,:) = mmr(:,:,id_CFC11) |
---|
| 440 | ! cfc12_inca(:,:) = mmr(:,:,id_CFC12) |
---|
| 441 | ! h2o_inca(:,:) = mmr(:,:,id_H2O) |
---|
| 442 | ! |
---|
| 443 | ! |
---|
| 444 | ! |
---|
| 445 | ! ! Ozone field POZON (kg/kg) |
---|
| 446 | ! ! Note: For POZON (used in radlwsw) simply use o3_inca from above (in mmr) |
---|
| 447 | ! |
---|
| 448 | ! ! 2D fields for CH4, N2O, CFC11 and CFC12 (in mmr) |
---|
| 449 | ! DO jj = 1 , PLEV |
---|
| 450 | ! jki=(jj-1)*ng1p1+1 |
---|
| 451 | ! jkf=jki+ng1 |
---|
| 452 | ! DO jk = jki, jkf |
---|
| 453 | ! ch4_inca_2d(:,jk) = ch4_inca(:,jj) |
---|
| 454 | ! n2o_inca_2d(:,jk) = n2o_inca(:,jj) |
---|
| 455 | ! cfc11_inca_2d(:,jk) = cfc11_inca(:,jj) |
---|
| 456 | ! cfc12_inca_2d(:,jk) = cfc12_inca(:,jj) |
---|
| 457 | ! ENDDO |
---|
| 458 | ! ENDDO |
---|
| 459 | |
---|
| 460 | END SUBROUTINE FIELD_PREP |
---|
| 461 | |
---|
| 462 | |
---|
| 463 | SUBROUTINE OZLIN_READ(calday) |
---|
| 464 | !---------------------------------------------------------------------- |
---|
| 465 | ! ... Read ozone linear coefficients |
---|
| 466 | ! Didier Hauglustaine, IPSL, 1999. |
---|
| 467 | !---------------------------------------------------------------------- |
---|
| 468 | |
---|
| 469 | USE O3LIN_COM |
---|
| 470 | USE INCA_DIM |
---|
| 471 | USE MOD_INCA_PARA |
---|
| 472 | USE PRINT_INCA |
---|
| 473 | |
---|
| 474 | IMPLICIT NONE |
---|
| 475 | |
---|
| 476 | INCLUDE 'netcdf.inc' |
---|
| 477 | |
---|
| 478 | !---------------------------------------------------------------------- |
---|
| 479 | ! ... Dummy args |
---|
| 480 | !---------------------------------------------------------------------- |
---|
| 481 | REAL, INTENT(IN) :: calday |
---|
| 482 | |
---|
| 483 | !---------------------------------------------------------------------- |
---|
| 484 | ! ... Local variables |
---|
| 485 | !---------------------------------------------------------------------- |
---|
| 486 | INTEGER :: iret |
---|
| 487 | INTEGER :: varid1,varid2,varid3,varid4,varid5 |
---|
| 488 | INTEGER :: varid6,varid7,varid8,varid9 |
---|
| 489 | INTEGER :: oldlotime, oldhitime |
---|
| 490 | INTEGER, DIMENSION(3) :: start, count |
---|
| 491 | REAL :: A1bd_glo(nbp_glo,nlevl,2) |
---|
| 492 | REAL :: A2bd_glo(nbp_glo,nlevl,2) |
---|
| 493 | REAL :: A3bd_glo(nbp_glo,nlevl,2) |
---|
| 494 | REAL :: A4bd_glo(nbp_glo,nlevl,2) |
---|
| 495 | REAL :: A5bd_glo(nbp_glo,nlevl,2) |
---|
| 496 | REAL :: A6bd_glo(nbp_glo,nlevl,2) |
---|
| 497 | REAL :: A7bd_glo(nbp_glo,nlevl,2) |
---|
| 498 | REAL :: A8bd_glo(nbp_glo,nlevl,2) |
---|
| 499 | REAL :: A9bd_glo(nbp_glo,nlevl,2) |
---|
| 500 | |
---|
| 501 | ! ... Check to see if model time still bounded by data set |
---|
| 502 | oldlotime = lotime |
---|
| 503 | oldhitime = hitime |
---|
| 504 | CALL FINDPLB (timel, ntimel, calday, lotime) |
---|
| 505 | IF ( cyclical ) THEN |
---|
| 506 | hitime = MOD(lotime,ntimel) + 1 |
---|
| 507 | ELSE |
---|
| 508 | hitime = lotime + 1 |
---|
| 509 | END IF |
---|
| 510 | |
---|
| 511 | ! ... Read new data if needed |
---|
| 512 | IF ( hitime /= oldhitime ) THEN |
---|
| 513 | loind = hiind |
---|
| 514 | hiind = MOD( loind, 2) + 1 |
---|
| 515 | |
---|
| 516 | !$OMP MASTER |
---|
| 517 | IF (is_mpi_root) THEN |
---|
| 518 | iret = nf_inq_varid(ncidl, 'A1', varid1) |
---|
| 519 | CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A1') |
---|
| 520 | iret = nf_inq_varid(ncidl, 'A2', varid2) |
---|
| 521 | CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A2') |
---|
| 522 | iret = nf_inq_varid(ncidl, 'A3', varid3) |
---|
| 523 | CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A3') |
---|
| 524 | iret = nf_inq_varid(ncidl, 'A4', varid4) |
---|
| 525 | CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A4') |
---|
| 526 | iret = nf_inq_varid(ncidl, 'A5', varid5) |
---|
| 527 | CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A5') |
---|
| 528 | iret = nf_inq_varid(ncidl, 'A6', varid6) |
---|
| 529 | CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A6') |
---|
| 530 | iret = nf_inq_varid(ncidl, 'A7', varid7) |
---|
| 531 | CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A7') |
---|
| 532 | iret = nf_inq_varid(ncidl, 'A8', varid8) |
---|
| 533 | CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A8') |
---|
| 534 | iret = nf_inq_varid(ncidl, 'A9', varid9) |
---|
| 535 | CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A9') |
---|
| 536 | |
---|
| 537 | count(1) = nbp_glo |
---|
| 538 | count(2) = nlevl |
---|
| 539 | count(3) = 1 |
---|
| 540 | start(1) = 1 |
---|
| 541 | start(2) = 1 |
---|
| 542 | |
---|
| 543 | WRITE(lunout,*) 'OZLIN_READ : read new data for time ', timel(hitime) |
---|
| 544 | start(3) = hitime |
---|
| 545 | iret = nf_get_vara_double(ncidl, varid1, start, count, A1bd_glo(1,1,hiind)) |
---|
| 546 | CALL check_err(iret, 'OCLIN_READ', 'problem to read A1bd for hiind') |
---|
| 547 | iret = nf_get_vara_double(ncidl, varid2, start, count, A2bd_glo(1,1,hiind)) |
---|
| 548 | CALL check_err(iret, 'OCLIN_READ', 'problem to read A2bd for hiind') |
---|
| 549 | iret = nf_get_vara_double(ncidl, varid3, start, count, A3bd_glo(1,1,hiind)) |
---|
| 550 | CALL check_err(iret, 'OCLIN_READ', 'problem to read A3bd for hiind') |
---|
| 551 | iret = nf_get_vara_double(ncidl, varid4, start, count, A4bd_glo(1,1,hiind)) |
---|
| 552 | CALL check_err(iret, 'OCLIN_READ', 'problem to read A4bd for hiind') |
---|
| 553 | iret = nf_get_vara_double(ncidl, varid5, start, count, A5bd_glo(1,1,hiind)) |
---|
| 554 | CALL check_err(iret, 'OCLIN_READ', 'problem to read A5bd for hiind') |
---|
| 555 | iret = nf_get_vara_double(ncidl, varid6, start, count, A6bd_glo(1,1,hiind)) |
---|
| 556 | CALL check_err(iret, 'OCLIN_READ', 'problem to read A6bd for hiind') |
---|
| 557 | iret = nf_get_vara_double(ncidl, varid7, start, count, A7bd_glo(1,1,hiind)) |
---|
| 558 | CALL check_err(iret, 'OCLIN_READ', 'problem to read A7bd for hiind') |
---|
| 559 | iret = nf_get_vara_double(ncidl, varid8, start, count, A8bd_glo(1,1,hiind)) |
---|
| 560 | CALL check_err(iret, 'OCLIN_READ', 'problem to read A8bd for hiind') |
---|
| 561 | iret = nf_get_vara_double(ncidl, varid9, start, count, A9bd_glo(1,1,hiind)) |
---|
| 562 | CALL check_err(iret, 'OCLIN_READ', 'problem to read A9bd for hiind') |
---|
| 563 | |
---|
| 564 | ENDIF |
---|
| 565 | !$OMP END MASTER |
---|
| 566 | |
---|
| 567 | CALL scatter(A1bd_glo(:,:,hiind),A1bd(:,:,hiind)) |
---|
| 568 | CALL scatter(A2bd_glo(:,:,hiind),A2bd(:,:,hiind)) |
---|
| 569 | CALL scatter(A3bd_glo(:,:,hiind),A3bd(:,:,hiind)) |
---|
| 570 | CALL scatter(A4bd_glo(:,:,hiind),A4bd(:,:,hiind)) |
---|
| 571 | CALL scatter(A5bd_glo(:,:,hiind),A5bd(:,:,hiind)) |
---|
| 572 | CALL scatter(A6bd_glo(:,:,hiind),A6bd(:,:,hiind)) |
---|
| 573 | CALL scatter(A7bd_glo(:,:,hiind),A7bd(:,:,hiind)) |
---|
| 574 | CALL scatter(A8bd_glo(:,:,hiind),A8bd(:,:,hiind)) |
---|
| 575 | CALL scatter(A9bd_glo(:,:,hiind),A9bd(:,:,hiind)) |
---|
| 576 | |
---|
| 577 | IF ( lotime /= oldhitime ) THEN |
---|
| 578 | !$OMP MASTER |
---|
| 579 | IF (is_mpi_root) THEN |
---|
| 580 | WRITE(lunout,*) 'OZLIN_READ : read new data for time ', timel(lotime) |
---|
| 581 | start(3) = lotime |
---|
| 582 | iret = nf_get_vara_double(ncidl, varid1, start, count, A1bd_glo(1,1,loind)) |
---|
| 583 | CALL check_err(iret, 'OZLIN_READ', 'problem to read A1bd for loind') |
---|
| 584 | iret = nf_get_vara_double(ncidl, varid2, start, count, A2bd_glo(1,1,loind)) |
---|
| 585 | CALL check_err(iret, 'OZLIN_READ', 'problem to read A2bd for loind') |
---|
| 586 | iret = nf_get_vara_double(ncidl, varid3, start, count, A3bd_glo(1,1,loind)) |
---|
| 587 | CALL check_err(iret, 'OZLIN_READ', 'problem to read A3bd for loind') |
---|
| 588 | iret = nf_get_vara_double(ncidl, varid4, start, count, A4bd_glo(1,1,loind)) |
---|
| 589 | CALL check_err(iret, 'OZLIN_READ', 'problem to read A4bd for loind') |
---|
| 590 | iret = nf_get_vara_double(ncidl, varid5, start, count, A5bd_glo(1,1,loind)) |
---|
| 591 | CALL check_err(iret, 'OZLIN_READ', 'problem to read A5bd for loind') |
---|
| 592 | iret = nf_get_vara_double(ncidl, varid6, start, count, A6bd_glo(1,1,loind)) |
---|
| 593 | CALL check_err(iret, 'OZLIN_READ', 'problem to read A6bd for loind') |
---|
| 594 | iret = nf_get_vara_double(ncidl, varid7, start, count, A7bd_glo(1,1,loind)) |
---|
| 595 | CALL check_err(iret, 'OZLIN_READ', 'problem to read A7bd for loind') |
---|
| 596 | iret = nf_get_vara_double(ncidl, varid8, start, count, A8bd_glo(1,1,loind)) |
---|
| 597 | CALL check_err(iret, 'OZLIN_READ', 'problem to read A8bd for loind') |
---|
| 598 | iret = nf_get_vara_double(ncidl, varid9, start, count, A9bd_glo(1,1,loind)) |
---|
| 599 | CALL check_err(iret, 'OZLIN_READ', 'problem to read A9bd for loind') |
---|
| 600 | |
---|
| 601 | ENDIF |
---|
| 602 | !$OMP END MASTER |
---|
| 603 | CALL scatter(A1bd_glo(:,:,loind),A1bd(:,:,loind)) |
---|
| 604 | CALL scatter(A2bd_glo(:,:,loind),A2bd(:,:,loind)) |
---|
| 605 | CALL scatter(A3bd_glo(:,:,loind),A3bd(:,:,loind)) |
---|
| 606 | CALL scatter(A4bd_glo(:,:,loind),A4bd(:,:,loind)) |
---|
| 607 | CALL scatter(A5bd_glo(:,:,loind),A5bd(:,:,loind)) |
---|
| 608 | CALL scatter(A6bd_glo(:,:,loind),A6bd(:,:,loind)) |
---|
| 609 | CALL scatter(A7bd_glo(:,:,loind),A7bd(:,:,loind)) |
---|
| 610 | CALL scatter(A8bd_glo(:,:,loind),A8bd(:,:,loind)) |
---|
| 611 | CALL scatter(A9bd_glo(:,:,loind),A9bd(:,:,loind)) |
---|
| 612 | |
---|
| 613 | END IF |
---|
| 614 | |
---|
| 615 | END IF |
---|
| 616 | |
---|
| 617 | END SUBROUTINE OZLIN_READ |
---|
| 618 | |
---|
| 619 | SUBROUTINE OZLIN_INTERP(calday, pmid) |
---|
| 620 | !---------------------------------------------------------------------- |
---|
| 621 | ! ... Interpolate ozone linear coefficients on model time and on pressure levels |
---|
| 622 | !R. Valorso |
---|
| 623 | !---------------------------------------------------------------------- |
---|
| 624 | |
---|
| 625 | USE O3LIN_COM |
---|
| 626 | USE INCA_DIM |
---|
| 627 | IMPLICIT NONE |
---|
| 628 | |
---|
| 629 | !---------------------------------------------------------------------- |
---|
| 630 | ! ... Dummy args |
---|
| 631 | !---------------------------------------------------------------------- |
---|
| 632 | REAL, INTENT(IN) :: calday |
---|
| 633 | REAL, INTENT(IN) :: pmid(PLON,PLEV) |
---|
| 634 | |
---|
| 635 | !---------------------------------------------------------------------- |
---|
| 636 | ! ... Local variables |
---|
| 637 | !---------------------------------------------------------------------- |
---|
| 638 | REAL :: dt, dt1, tint |
---|
| 639 | REAL :: ploc(klonl,nlevl) |
---|
| 640 | INTEGER :: lonlev |
---|
| 641 | INTEGER :: i |
---|
| 642 | INTEGER :: index(PLON,PLEV) |
---|
| 643 | |
---|
| 644 | index = -9999 |
---|
| 645 | lonlev = klonl * nlevl |
---|
| 646 | |
---|
| 647 | DO i = 1, klonl |
---|
| 648 | ploc(i,:) = presl(:) |
---|
| 649 | END DO |
---|
| 650 | |
---|
| 651 | !---------------------------------------------------------------------- |
---|
| 652 | ! Note: 365 day year inconsistent with LMDz (360 days) !!! |
---|
| 653 | !---------------------------------------------------------------------- |
---|
| 654 | ! ... First interpolate linearly on current model time |
---|
| 655 | |
---|
| 656 | IF ( timel(hitime) < timel(lotime) ) THEN |
---|
| 657 | dt = 365. - timel(lotime) + timel(hitime) |
---|
| 658 | IF (calday <= timel(hitime) ) THEN |
---|
| 659 | dt1 = 365. - timel(lotime) + calday |
---|
| 660 | ELSE |
---|
| 661 | dt1 = calday - timel(lotime) |
---|
| 662 | END IF |
---|
| 663 | ELSE |
---|
| 664 | dt = timel(hitime) - timel(lotime) |
---|
| 665 | dt1 = calday - timel(lotime) |
---|
| 666 | END IF |
---|
| 667 | |
---|
| 668 | tint = dt1/dt |
---|
| 669 | |
---|
| 670 | CALL linintp (lonlev, & |
---|
| 671 | 0., & |
---|
| 672 | 1., & |
---|
| 673 | tint, & |
---|
| 674 | A1bd(1,1,loind), & |
---|
| 675 | A1bd(1,1,hiind), & |
---|
| 676 | A1cr) |
---|
| 677 | |
---|
| 678 | ! ... Put on model pressure levels |
---|
| 679 | |
---|
| 680 | CALL pinterp (PLON, & |
---|
| 681 | A1cr, & |
---|
| 682 | ploc, & |
---|
| 683 | nlevl, & |
---|
| 684 | A1, & |
---|
| 685 | pmid, & |
---|
| 686 | PLEV, & |
---|
| 687 | index, & |
---|
| 688 | .false. ) |
---|
| 689 | |
---|
| 690 | ! ... |
---|
| 691 | |
---|
| 692 | CALL linintp (lonlev, & |
---|
| 693 | 0., & |
---|
| 694 | 1., & |
---|
| 695 | tint, & |
---|
| 696 | A2bd(1,1,loind), & |
---|
| 697 | A2bd(1,1,hiind), & |
---|
| 698 | A2cr) |
---|
| 699 | |
---|
| 700 | ! ... Put on model pressure levels |
---|
| 701 | |
---|
| 702 | index = -9999 |
---|
| 703 | CALL pinterp (PLON, & |
---|
| 704 | A2cr, & |
---|
| 705 | ploc, & |
---|
| 706 | nlevl, & |
---|
| 707 | A2, & |
---|
| 708 | pmid, & |
---|
| 709 | PLEV, & |
---|
| 710 | index, & |
---|
| 711 | .false. ) |
---|
| 712 | ! ... |
---|
| 713 | |
---|
| 714 | CALL linintp (lonlev, & |
---|
| 715 | 0., & |
---|
| 716 | 1., & |
---|
| 717 | tint, & |
---|
| 718 | A3bd(1,1,loind), & |
---|
| 719 | A3bd(1,1,hiind), & |
---|
| 720 | A3cr) |
---|
| 721 | |
---|
| 722 | ! ... Put on model pressure levels |
---|
| 723 | |
---|
| 724 | index = -9999 |
---|
| 725 | CALL pinterp (PLON, & |
---|
| 726 | A3cr, & |
---|
| 727 | ploc, & |
---|
| 728 | nlevl, & |
---|
| 729 | A3, & |
---|
| 730 | pmid, & |
---|
| 731 | PLEV, & |
---|
| 732 | index, & |
---|
| 733 | .false. ) |
---|
| 734 | ! ... |
---|
| 735 | |
---|
| 736 | CALL linintp (lonlev, & |
---|
| 737 | 0., & |
---|
| 738 | 1., & |
---|
| 739 | tint, & |
---|
| 740 | A4bd(1,1,loind), & |
---|
| 741 | A4bd(1,1,hiind), & |
---|
| 742 | A4cr) |
---|
| 743 | |
---|
| 744 | ! ... Put on model pressure levels |
---|
| 745 | |
---|
| 746 | index = -9999 |
---|
| 747 | CALL pinterp (PLON, & |
---|
| 748 | A4cr, & |
---|
| 749 | ploc, & |
---|
| 750 | nlevl, & |
---|
| 751 | A4, & |
---|
| 752 | pmid, & |
---|
| 753 | PLEV, & |
---|
| 754 | index, & |
---|
| 755 | .false. ) |
---|
| 756 | ! ... |
---|
| 757 | |
---|
| 758 | CALL linintp (lonlev, & |
---|
| 759 | 0., & |
---|
| 760 | 1., & |
---|
| 761 | tint, & |
---|
| 762 | A5bd(1,1,loind), & |
---|
| 763 | A5bd(1,1,hiind), & |
---|
| 764 | A5cr) |
---|
| 765 | |
---|
| 766 | ! ... Put on model pressure levels |
---|
| 767 | |
---|
| 768 | index = -9999 |
---|
| 769 | CALL pinterp (PLON, & |
---|
| 770 | A5cr, & |
---|
| 771 | ploc, & |
---|
| 772 | nlevl, & |
---|
| 773 | A5, & |
---|
| 774 | pmid, & |
---|
| 775 | PLEV, & |
---|
| 776 | index, & |
---|
| 777 | .false. ) |
---|
| 778 | ! ... |
---|
| 779 | |
---|
| 780 | CALL linintp (lonlev, & |
---|
| 781 | 0., & |
---|
| 782 | 1., & |
---|
| 783 | tint, & |
---|
| 784 | A6bd(1,1,loind), & |
---|
| 785 | A6bd(1,1,hiind), & |
---|
| 786 | A6cr) |
---|
| 787 | |
---|
| 788 | ! ... Put on model pressure levels |
---|
| 789 | |
---|
| 790 | index = -9999 |
---|
| 791 | CALL pinterp (PLON, & |
---|
| 792 | A6cr, & |
---|
| 793 | ploc, & |
---|
| 794 | nlevl, & |
---|
| 795 | A6, & |
---|
| 796 | pmid, & |
---|
| 797 | PLEV, & |
---|
| 798 | index, & |
---|
| 799 | .false. ) |
---|
| 800 | ! ... |
---|
| 801 | |
---|
| 802 | CALL linintp (lonlev, & |
---|
| 803 | 0., & |
---|
| 804 | 1., & |
---|
| 805 | tint, & |
---|
| 806 | A7bd(1,1,loind), & |
---|
| 807 | A7bd(1,1,hiind), & |
---|
| 808 | A7cr) |
---|
| 809 | |
---|
| 810 | ! ... Put on model pressure levels |
---|
| 811 | |
---|
| 812 | index = -9999 |
---|
| 813 | CALL pinterp (PLON, & |
---|
| 814 | A7cr, & |
---|
| 815 | ploc, & |
---|
| 816 | nlevl, & |
---|
| 817 | A7, & |
---|
| 818 | pmid, & |
---|
| 819 | PLEV, & |
---|
| 820 | index, & |
---|
| 821 | .false. ) |
---|
| 822 | ! ... |
---|
| 823 | |
---|
| 824 | CALL linintp (lonlev, & |
---|
| 825 | 0., & |
---|
| 826 | 1., & |
---|
| 827 | tint, & |
---|
| 828 | A8bd(1,1,loind), & |
---|
| 829 | A8bd(1,1,hiind), & |
---|
| 830 | A8cr) |
---|
| 831 | |
---|
| 832 | ! ... Put on model pressure levels |
---|
| 833 | |
---|
| 834 | index = -9999 |
---|
| 835 | CALL pinterp (PLON, & |
---|
| 836 | A8cr, & |
---|
| 837 | ploc, & |
---|
| 838 | nlevl, & |
---|
| 839 | A8, & |
---|
| 840 | pmid, & |
---|
| 841 | PLEV, & |
---|
| 842 | index, & |
---|
| 843 | .false. ) |
---|
| 844 | ! ... |
---|
| 845 | |
---|
| 846 | CALL linintp (lonlev, & |
---|
| 847 | 0., & |
---|
| 848 | 1., & |
---|
| 849 | tint, & |
---|
| 850 | A9bd(1,1,loind), & |
---|
| 851 | A9bd(1,1,hiind), & |
---|
| 852 | A9cr) |
---|
| 853 | |
---|
| 854 | ! ... Put on model pressure levels |
---|
| 855 | |
---|
| 856 | index = -9999 |
---|
| 857 | CALL pinterp (PLON, & |
---|
| 858 | A9cr, & |
---|
| 859 | ploc, & |
---|
| 860 | nlevl, & |
---|
| 861 | A9, & |
---|
| 862 | pmid, & |
---|
| 863 | PLEV, & |
---|
| 864 | index, & |
---|
| 865 | .false. ) |
---|
| 866 | |
---|
| 867 | END SUBROUTINE OZLIN_INTERP |
---|
| 868 | |
---|
| 869 | |
---|
| 870 | #ifdef STRAT |
---|
| 871 | SUBROUTINE SAD_READ(calday) |
---|
| 872 | !---------------------------------------------------------------------- |
---|
| 873 | ! ... Read sulfate data from CMIP5 |
---|
| 874 | ! R. Valorso |
---|
| 875 | !---------------------------------------------------------------------- |
---|
| 876 | |
---|
| 877 | USE SAD_COM |
---|
| 878 | USE INCA_DIM |
---|
| 879 | USE MOD_INCA_PARA |
---|
| 880 | USE PRINT_INCA |
---|
| 881 | |
---|
| 882 | IMPLICIT NONE |
---|
| 883 | |
---|
| 884 | INCLUDE 'netcdf.inc' |
---|
| 885 | |
---|
| 886 | !---------------------------------------------------------------------- |
---|
| 887 | ! ... Dummy args |
---|
| 888 | !---------------------------------------------------------------------- |
---|
| 889 | REAL, INTENT(IN) :: calday |
---|
| 890 | |
---|
| 891 | !---------------------------------------------------------------------- |
---|
| 892 | ! ... Local variables |
---|
| 893 | !---------------------------------------------------------------------- |
---|
| 894 | INTEGER :: iret |
---|
| 895 | INTEGER :: varid1,varid2,varid3,varid4,varid5 |
---|
| 896 | INTEGER :: oldlotime, oldhitime |
---|
| 897 | INTEGER, DIMENSION(3) :: start, count |
---|
| 898 | REAL :: PRES_bd_glo(nbp_glo,nlevs,2) |
---|
| 899 | REAL :: SAD_bd_glo(nbp_glo,nlevs,2) |
---|
| 900 | REAL :: MASS_bd_glo(nbp_glo,nlevs,2) |
---|
| 901 | REAL :: VOL_bd_glo(nbp_glo,nlevs,2) |
---|
| 902 | REAL :: RMEAN_bd_glo(nbp_glo,nlevs,2) |
---|
| 903 | |
---|
| 904 | ! ... Check to see if model time still bounded by data set |
---|
| 905 | oldlotime = lotime |
---|
| 906 | oldhitime = hitime |
---|
| 907 | CALL FINDPLB (times, ntimes, calday, lotime) |
---|
| 908 | IF ( cyclical ) THEN |
---|
| 909 | hitime = MOD(lotime,ntimes) + 1 |
---|
| 910 | ELSE |
---|
| 911 | hitime = lotime + 1 |
---|
| 912 | END IF |
---|
| 913 | |
---|
| 914 | ! ... Read new data if needed |
---|
| 915 | IF ( hitime /= oldhitime ) THEN |
---|
| 916 | loind = hiind |
---|
| 917 | hiind = MOD( loind, 2) + 1 |
---|
| 918 | |
---|
| 919 | !$OMP MASTER |
---|
| 920 | IF (is_mpi_root) THEN |
---|
| 921 | iret = nf_inq_varid(ncids, 'PRES', varid1) |
---|
| 922 | CALL check_err(iret, 'SAF_READ', 'problem to get id for pres') |
---|
| 923 | iret = nf_inq_varid(ncids, 'SAD', varid2) |
---|
| 924 | CALL check_err(iret, 'SAD_READ', 'problem to get id for sad') |
---|
| 925 | iret = nf_inq_varid(ncids, 'MASS', varid3) |
---|
| 926 | CALL check_err(iret, 'SAD_READ', 'problem to get id for mass') |
---|
| 927 | iret = nf_inq_varid(ncids, 'VOL', varid4) |
---|
| 928 | CALL check_err(iret, 'SAD_READ', 'problem to get id for vol') |
---|
| 929 | iret = nf_inq_varid(ncids, 'RMEAN', varid5) |
---|
| 930 | CALL check_err(iret, 'SAD_READ', 'problem to get id for rmean') |
---|
| 931 | |
---|
| 932 | count(1) = nbp_glo |
---|
| 933 | count(2) = nlevs |
---|
| 934 | count(3) = 1 |
---|
| 935 | start(1) = 1 |
---|
| 936 | start(2) = 1 |
---|
| 937 | |
---|
| 938 | start(3) = hitime |
---|
| 939 | iret = nf_get_vara_double(ncids, varid1, start, count, PRES_bd_glo(1,1,hiind)) |
---|
| 940 | CALL check_err(iret, 'SAD_READ', 'problem to read PRES_bd hiind') |
---|
| 941 | iret = nf_get_vara_double(ncids, varid2, start, count, SAD_bd_glo(1,1,hiind)) |
---|
| 942 | CALL check_err(iret, 'SAD_READ', 'problem to read SAD_bd hiind') |
---|
| 943 | iret = nf_get_vara_double(ncids, varid3, start, count, MASS_bd_glo(1,1,hiind)) |
---|
| 944 | CALL check_err(iret, 'SAD_READ', 'problem to read MASS_bd hiind') |
---|
| 945 | iret = nf_get_vara_double(ncids, varid4, start, count, VOL_bd_glo(1,1,hiind)) |
---|
| 946 | CALL check_err(iret, 'SAD_READ', 'problem to read VOL_bd hiind') |
---|
| 947 | iret = nf_get_vara_double(ncids, varid5, start, count, RMEAN_bd_glo(1,1,hiind)) |
---|
| 948 | CALL check_err(iret, 'SAD_READ', 'problem to read RMEAN_bd hiind') |
---|
| 949 | |
---|
| 950 | ENDIF |
---|
| 951 | !$OMP END MASTER |
---|
| 952 | CALL scatter(PRES_bd_glo(:,:,hiind),PRES_bd(:,:,hiind)) |
---|
| 953 | CALL scatter(SAD_bd_glo(:,:,hiind),SAD_bd(:,:,hiind)) |
---|
| 954 | CALL scatter(MASS_bd_glo(:,:,hiind),MASS_bd(:,:,hiind)) |
---|
| 955 | CALL scatter(VOL_bd_glo(:,:,hiind),VOL_bd(:,:,hiind)) |
---|
| 956 | CALL scatter(RMEAN_bd_glo(:,:,hiind),RMEAN_bd(:,:,hiind)) |
---|
| 957 | |
---|
| 958 | IF ( lotime /= oldhitime ) THEN |
---|
| 959 | !$OMP MASTER |
---|
| 960 | IF (is_mpi_root) THEN |
---|
| 961 | start(3) = lotime |
---|
| 962 | iret = nf_get_vara_double(ncids, varid1, start, count, PRES_bd_glo(1,1,loind)) |
---|
| 963 | CALL check_err(iret, 'SAD_READ', 'problem to read PRES_bd loind') |
---|
| 964 | iret = nf_get_vara_double(ncids, varid2, start, count, SAD_bd_glo(1,1,loind)) |
---|
| 965 | CALL check_err(iret, 'SAD_READ', 'problem to read SAD_bd loind') |
---|
| 966 | iret = nf_get_vara_double(ncids, varid3, start, count, MASS_bd_glo(1,1,loind)) |
---|
| 967 | CALL check_err(iret, 'SAD_READ', 'problem to read MASS_bd loind') |
---|
| 968 | iret = nf_get_vara_double(ncids, varid4, start, count, VOL_bd_glo(1,1,loind)) |
---|
| 969 | CALL check_err(iret, 'SAD_READ', 'problem to read VOL_bd loind') |
---|
| 970 | iret = nf_get_vara_double(ncids, varid5, start, count, RMEAN_bd_glo(1,1,loind)) |
---|
| 971 | CALL check_err(iret, 'SAD_READ', 'problem to read RMEAN_bd loind') |
---|
| 972 | |
---|
| 973 | ENDIF |
---|
| 974 | !$OMP END MASTER |
---|
| 975 | CALL scatter(PRES_bd_glo(:,:,loind),PRES_bd(:,:,loind)) |
---|
| 976 | CALL scatter(SAD_bd_glo(:,:,loind),SAD_bd(:,:,loind)) |
---|
| 977 | CALL scatter(MASS_bd_glo(:,:,loind),MASS_bd(:,:,loind)) |
---|
| 978 | CALL scatter(VOL_bd_glo(:,:,loind),VOL_bd(:,:,loind)) |
---|
| 979 | CALL scatter(RMEAN_bd_glo(:,:,loind),RMEAN_bd(:,:,loind)) |
---|
| 980 | |
---|
| 981 | END IF |
---|
| 982 | |
---|
| 983 | END IF |
---|
| 984 | |
---|
| 985 | END SUBROUTINE SAD_READ |
---|
| 986 | |
---|
| 987 | SUBROUTINE SAD_INTERP(calday, pmid) |
---|
| 988 | !---------------------------------------------------------------------- |
---|
| 989 | ! ... Interpolate sulfate data from CMIP 5 on model time and on pressure levels |
---|
| 990 | ! R. Valorso |
---|
| 991 | !---------------------------------------------------------------------- |
---|
| 992 | |
---|
| 993 | USE SAD_COM |
---|
| 994 | USE INCA_DIM |
---|
| 995 | IMPLICIT NONE |
---|
| 996 | |
---|
| 997 | !---------------------------------------------------------------------- |
---|
| 998 | ! ... Dummy args |
---|
| 999 | !---------------------------------------------------------------------- |
---|
| 1000 | REAL, INTENT(IN) :: calday |
---|
| 1001 | REAL, INTENT(IN) :: pmid(PLON,PLEV) |
---|
| 1002 | |
---|
| 1003 | !---------------------------------------------------------------------- |
---|
| 1004 | ! ... Local variables |
---|
| 1005 | !---------------------------------------------------------------------- |
---|
| 1006 | REAL :: dt, dt1, tint |
---|
| 1007 | REAL :: ploc(klons,nlevs) |
---|
| 1008 | INTEGER :: lonlev |
---|
| 1009 | INTEGER :: i |
---|
| 1010 | INTEGER :: index(PLON,PLEV) |
---|
| 1011 | |
---|
| 1012 | index = -9999 |
---|
| 1013 | lonlev = klons * nlevs |
---|
| 1014 | |
---|
| 1015 | !---------------------------------------------------------------------- |
---|
| 1016 | ! ... First interpolate linearly on current model time |
---|
| 1017 | |
---|
| 1018 | IF ( times(hitime) < times(lotime) ) THEN |
---|
| 1019 | dt = 365. - times(lotime) + times(hitime) |
---|
| 1020 | IF (calday <= times(hitime) ) THEN |
---|
| 1021 | dt1 = 365. - times(lotime) + calday |
---|
| 1022 | ELSE |
---|
| 1023 | dt1 = calday - times(lotime) |
---|
| 1024 | END IF |
---|
| 1025 | ELSE |
---|
| 1026 | dt = times(hitime) - times(lotime) |
---|
| 1027 | dt1 = calday - times(lotime) |
---|
| 1028 | END IF |
---|
| 1029 | |
---|
| 1030 | tint = dt1/dt |
---|
| 1031 | |
---|
| 1032 | CALL linintp (lonlev, & |
---|
| 1033 | 0., & |
---|
| 1034 | 1., & |
---|
| 1035 | tint, & |
---|
| 1036 | PRES_bd(1,1,loind), & |
---|
| 1037 | PRES_bd(1,1,hiind), & |
---|
| 1038 | PRES_cr) |
---|
| 1039 | |
---|
| 1040 | DO i = 1, klons |
---|
| 1041 | ploc(i,:) = PRES_cr(i,:) |
---|
| 1042 | ENDDO |
---|
| 1043 | |
---|
| 1044 | ! ... |
---|
| 1045 | |
---|
| 1046 | CALL linintp (lonlev, & |
---|
| 1047 | 0., & |
---|
| 1048 | 1., & |
---|
| 1049 | tint, & |
---|
| 1050 | SAD_bd(1,1,loind), & |
---|
| 1051 | SAD_bd(1,1,hiind), & |
---|
| 1052 | SAD_cr) |
---|
| 1053 | |
---|
| 1054 | ! ... Put on model pressure levels |
---|
| 1055 | |
---|
| 1056 | index = -9999 |
---|
| 1057 | CALL pinterp (PLON, & |
---|
| 1058 | SAD_cr, & |
---|
| 1059 | ploc, & |
---|
| 1060 | nlevs, & |
---|
| 1061 | SAD, & |
---|
| 1062 | pmid, & |
---|
| 1063 | PLEV, & |
---|
| 1064 | index, & |
---|
| 1065 | .false. ) |
---|
| 1066 | ! ... |
---|
| 1067 | |
---|
| 1068 | CALL linintp (lonlev, & |
---|
| 1069 | 0., & |
---|
| 1070 | 1., & |
---|
| 1071 | tint, & |
---|
| 1072 | MASS_bd(1,1,loind), & |
---|
| 1073 | MASS_bd(1,1,hiind), & |
---|
| 1074 | MASS_cr) |
---|
| 1075 | |
---|
| 1076 | ! ... Put on model pressure levels |
---|
| 1077 | |
---|
| 1078 | index = -9999 |
---|
| 1079 | CALL pinterp (PLON, & |
---|
| 1080 | MASS_cr, & |
---|
| 1081 | ploc, & |
---|
| 1082 | nlevs, & |
---|
| 1083 | MASS, & |
---|
| 1084 | pmid, & |
---|
| 1085 | PLEV, & |
---|
| 1086 | index, & |
---|
| 1087 | .false. ) |
---|
| 1088 | |
---|
| 1089 | ! ... |
---|
| 1090 | |
---|
| 1091 | CALL linintp (lonlev, & |
---|
| 1092 | 0., & |
---|
| 1093 | 1., & |
---|
| 1094 | tint, & |
---|
| 1095 | VOL_bd(1,1,loind), & |
---|
| 1096 | VOL_bd(1,1,hiind), & |
---|
| 1097 | VOL_cr) |
---|
| 1098 | |
---|
| 1099 | ! ... Put on model pressure levels |
---|
| 1100 | |
---|
| 1101 | index = -9999 |
---|
| 1102 | CALL pinterp (PLON, & |
---|
| 1103 | VOL_cr, & |
---|
| 1104 | ploc, & |
---|
| 1105 | nlevs, & |
---|
| 1106 | VOL, & |
---|
| 1107 | pmid, & |
---|
| 1108 | PLEV, & |
---|
| 1109 | index, & |
---|
| 1110 | .false. ) |
---|
| 1111 | ! ... |
---|
| 1112 | |
---|
| 1113 | CALL linintp (lonlev, & |
---|
| 1114 | 0., & |
---|
| 1115 | 1., & |
---|
| 1116 | tint, & |
---|
| 1117 | RMEAN_bd(1,1,loind), & |
---|
| 1118 | RMEAN_bd(1,1,hiind), & |
---|
| 1119 | RMEAN_cr) |
---|
| 1120 | |
---|
| 1121 | ! ... Put on model pressure levels |
---|
| 1122 | |
---|
| 1123 | index = -9999 |
---|
| 1124 | CALL pinterp (PLON, & |
---|
| 1125 | RMEAN_cr, & |
---|
| 1126 | ploc, & |
---|
| 1127 | nlevs, & |
---|
| 1128 | RMEAN, & |
---|
| 1129 | pmid, & |
---|
| 1130 | PLEV, & |
---|
| 1131 | index, & |
---|
| 1132 | .false. ) |
---|
| 1133 | |
---|
| 1134 | |
---|
| 1135 | END SUBROUTINE SAD_INTERP |
---|
| 1136 | #endif |
---|
| 1137 | #endif |
---|
| 1138 | #if !defined(GES) && !defined(DUSS) |
---|
| 1139 | SUBROUTINE BOUNDSPC (dtinv, mmr, pmid, temp) |
---|
| 1140 | !---------------------------------------------------------------------- |
---|
| 1141 | ! ... Impose boundary conditions for chemical tracers |
---|
| 1142 | ! Didier Hauglustaine, IPSL, 1999. |
---|
| 1143 | !---------------------------------------------------------------------- |
---|
| 1144 | |
---|
| 1145 | USE O3CLIM_COM |
---|
| 1146 | USE SPECIES_NAMES |
---|
| 1147 | USE OXYDANT_COM |
---|
| 1148 | USE AIRPLANE_SRC, ONLY : itrop |
---|
| 1149 | USE CHEM_MODS, ONLY : adv_mass, invariants |
---|
| 1150 | USE INCA_DIM |
---|
| 1151 | #ifdef STRAT |
---|
| 1152 | USE LGLIVED_MOD |
---|
| 1153 | #endif |
---|
| 1154 | USE PARAM_CHEM |
---|
| 1155 | USE RATE_INDEX_MOD |
---|
| 1156 | |
---|
| 1157 | IMPLICIT NONE |
---|
| 1158 | |
---|
| 1159 | !---------------------------------------------------------------------- |
---|
| 1160 | ! ... Dummy args |
---|
| 1161 | !---------------------------------------------------------------------- |
---|
| 1162 | REAL, INTENT(IN) :: dtinv |
---|
| 1163 | REAL, INTENT(IN) :: pmid(PLON,PLEV) |
---|
| 1164 | REAL, INTENT(IN) :: temp(PLON,PLEV) |
---|
| 1165 | REAL, INTENT(INOUT) :: mmr(PLON,PLEV,PCNST) |
---|
| 1166 | |
---|
| 1167 | !---------------------------------------------------------------------- |
---|
| 1168 | ! ... Local variables |
---|
| 1169 | !---------------------------------------------------------------------- |
---|
| 1170 | INTEGER :: i, k |
---|
| 1171 | INTEGER :: ktrop, krelax, k380 |
---|
| 1172 | REAL :: taurelax, facrelax |
---|
| 1173 | REAL :: delt |
---|
| 1174 | REAL, PARAMETER :: dry_mass = 28.966 |
---|
| 1175 | |
---|
| 1176 | REAL :: theta(PLON,PLEV) |
---|
| 1177 | |
---|
| 1178 | delt = 1./dtinv |
---|
| 1179 | taurelax = 86400.*10. ! (10 days) |
---|
| 1180 | facrelax = 1. - EXP( -delt / taurelax ) |
---|
| 1181 | theta(:,:) = temp(:,:)*(1.e5/pmid(:,:))**0.286 |
---|
| 1182 | |
---|
| 1183 | DO i = 1, PLON |
---|
| 1184 | |
---|
| 1185 | ktrop = itrop(i) |
---|
| 1186 | |
---|
| 1187 | DO k = ktrop, PLEV |
---|
| 1188 | IF (theta(i,k) >= 380.) THEN |
---|
| 1189 | k380 = k |
---|
| 1190 | EXIT |
---|
| 1191 | ENDIF |
---|
| 1192 | ENDDO |
---|
| 1193 | |
---|
| 1194 | ! krelax = MIN(ktrop+2,PLEV) ! 2 levels above tropopause |
---|
| 1195 | krelax = k380 ! 380K theta surface |
---|
| 1196 | |
---|
| 1197 | |
---|
| 1198 | ! ... Impose NOy to climatologies |
---|
| 1199 | ! mmr(i,krelax:PLEV,id_NO) = nomzrt(i,krelax:PLEV) * adv_mass(id_NO)/dry_mass |
---|
| 1200 | #ifdef AERONLY |
---|
| 1201 | mmr(i,krelax:PLEV,id_NO2) = invariants(i,krelax:PLEV, inv_NO2INV)/invariants(i,krelax:PLEV, inv_M) * adv_mass(id_NO2)/dry_mass |
---|
| 1202 | mmr(i,krelax:PLEV,id_HNO3) = invariants(i,krelax:PLEV, inv_HNO3INV)/invariants(i,krelax:PLEV, inv_M) * adv_mass(id_HNO3)/dry_mass |
---|
| 1203 | #else |
---|
| 1204 | ! ... Relax O3 to climatologies |
---|
| 1205 | if (trim(flag_o3) .eq. "o3clim") then |
---|
| 1206 | mmr(i,krelax:PLEV,id_O3) = o3clim(i,krelax:PLEV) |
---|
| 1207 | elseif (trim(flag_o3) .eq. "o3lin") then |
---|
| 1208 | mmr(i,krelax:PLEV,id_O3) = mmr(i,ktrop:PLEV,id_O3L) |
---|
| 1209 | endif |
---|
| 1210 | |
---|
| 1211 | ! ... Impose inert O3 to O3 above tropopause |
---|
| 1212 | mmr(i,ktrop:PLEV,id_O3I) = mmr(i,ktrop:PLEV,id_O3) |
---|
| 1213 | |
---|
| 1214 | ! ... Impose stratospheric O3 to O3 above tropopause |
---|
| 1215 | mmr(i,ktrop:PLEV,id_O3S) = mmr(i,ktrop:PLEV,id_O3) |
---|
| 1216 | |
---|
| 1217 | #endif |
---|
| 1218 | ! ... Impose CH4 |
---|
| 1219 | ! mmr(i,:,id_CH4) = 1760E-9 * 16./dry_mass ! forcage CH4 pour IPCC2000 |
---|
| 1220 | ! mmr(i,:,id_CH4) = 791.6E-9 * 16./dry_mass ! forcage CH4 pour IPCC2000 |
---|
| 1221 | |
---|
| 1222 | #ifdef STRAT |
---|
| 1223 | ! Stratospheric Version |
---|
| 1224 | ! Impose long_lived species at the surface CMIP5 |
---|
| 1225 | |
---|
| 1226 | mmr(i,1,id_CH4) = conc_cfc(1) * adv_mass(id_CH4) / dry_mass |
---|
| 1227 | mmr(i,1,id_N2O) = conc_cfc(2) * adv_mass(id_N2O) / dry_mass |
---|
| 1228 | |
---|
| 1229 | mmr(i,1,id_CFC11) = conc_cfc(3) * adv_mass(id_CFC11) / dry_mass |
---|
| 1230 | mmr(i,1,id_CFC12) = conc_cfc(4) * adv_mass(id_CFC12) / dry_mass |
---|
| 1231 | mmr(i,1,id_CFC113) = conc_cfc(5) * adv_mass(id_CFC113) / dry_mass |
---|
| 1232 | mmr(i,1,id_HCFC22) = conc_cfc(6) * adv_mass(id_HCFC22) / dry_mass |
---|
| 1233 | mmr(i,1,id_HFC134a) = conc_cfc(7) * adv_mass(id_HFC134a) / dry_mass |
---|
| 1234 | mmr(i,1,id_HCFC141) = conc_cfc(8) * adv_mass(id_HCFC141) / dry_mass |
---|
| 1235 | mmr(i,1,id_HCFC142) = conc_cfc(9) * adv_mass(id_HCFC142) / dry_mass |
---|
| 1236 | mmr(i,1,id_Ha1301) = conc_cfc(10) * adv_mass(id_Ha1301) / dry_mass |
---|
| 1237 | mmr(i,1,id_Ha1211) = conc_cfc(11) * adv_mass(id_Ha1211) / dry_mass |
---|
| 1238 | mmr(i,1,id_MCF) = conc_cfc(12) * adv_mass(id_MCF) / dry_mass |
---|
| 1239 | mmr(i,1,id_CCl4) = conc_cfc(13) * adv_mass(id_CCl4) / dry_mass |
---|
| 1240 | mmr(i,1,id_CH3Cl) = conc_cfc(14) * adv_mass(id_CH3Cl) / dry_mass |
---|
| 1241 | mmr(i,1,id_CH3Br) = conc_cfc(15) * adv_mass(id_CH3Br) / dry_mass |
---|
| 1242 | #endif |
---|
| 1243 | |
---|
| 1244 | END DO |
---|
| 1245 | |
---|
| 1246 | END SUBROUTINE BOUNDSPC |
---|
| 1247 | #endif |
---|
| 1248 | SUBROUTINE FDTROPOPAUSE (nx, ny, nz, pmid, temp) |
---|
| 1249 | !**************************************************************************** |
---|
| 1250 | ! Purpose |
---|
| 1251 | ! ------- |
---|
| 1252 | ! |
---|
| 1253 | ! Calculates 2D fields of thermal tropopause pressures and tropopause |
---|
| 1254 | ! level indices for given 3D fields of temperature and pressure. |
---|
| 1255 | ! |
---|
| 1256 | ! Methods |
---|
| 1257 | ! ------- |
---|
| 1258 | ! - Subroutine stattrop calculates the thermal tropopause pressure (Pa) |
---|
| 1259 | ! for a 1D column (sounding) of pressures (Pa) and temperatures (K). |
---|
| 1260 | ! |
---|
| 1261 | ! - The program tropo_test also corrects for non-sense occuring at high |
---|
| 1262 | ! latitudes > 60N (60S). |
---|
| 1263 | ! Theese problems occur at high latitudes on the winter hemisphere because |
---|
| 1264 | ! of the cold stratosphere and the consequently weak troposphere- |
---|
| 1265 | ! stratosphere transition of the thermal lapse rate -dT/dz which is |
---|
| 1266 | ! used to determine the thermal tropopause. |
---|
| 1267 | ! The problem is much less severe on the winter hemisphere. |
---|
| 1268 | ! |
---|
| 1269 | ! Here a fix is implemented applying a factor which scales the |
---|
| 1270 | ! tropopause pressures at latitudes LAT > 60N (60S) to the |
---|
| 1271 | ! mean tropopause pressure at 60N (60S), whenever the mean |
---|
| 1272 | ! tropopause pressure at LAT is lower than that at 60N (60S). |
---|
| 1273 | ! This idea was formulated in the ECHAM4/CHEM routine xttphwmo by |
---|
| 1274 | ! D. Nodorp, Ch. Land, B. Steil and R. Hein. |
---|
| 1275 | ! |
---|
| 1276 | ! Programmed by Dominik Brunner, ETHZ, Switzerland V.01, 10 Aug 2000 |
---|
| 1277 | ! Modified by Didier Hauglustaine, IPSL, for LMDZ/INCA, Oct 2000. |
---|
| 1278 | ! |
---|
| 1279 | !***************************************************************************** |
---|
| 1280 | |
---|
| 1281 | USE AIRPLANE_SRC, ONLY : itrop, ptrop |
---|
| 1282 | USE INCA_DIM |
---|
| 1283 | USE MOD_INCA_PARA, ONLY : ij_begin,ij_end, & |
---|
| 1284 | nbp_para_begin,nbp_para_end,mpi_rank, & |
---|
| 1285 | plon_omp_begin, plon_omp_end |
---|
| 1286 | USE MOD_GRID_INCA |
---|
| 1287 | USE MOD_INCA_OMP_DATA |
---|
| 1288 | |
---|
| 1289 | IMPLICIT NONE |
---|
| 1290 | !---------------------------------------------------------------------- |
---|
| 1291 | ! ... Dummy args |
---|
| 1292 | !---------------------------------------------------------------------- |
---|
| 1293 | REAL, INTENT(IN) :: pmid(PLON_GLO,PLEV) |
---|
| 1294 | REAL, INTENT(IN) :: temp(PLON_GLO,PLEV) |
---|
| 1295 | INTEGER, INTENT(IN) :: nx, ny, nz |
---|
| 1296 | |
---|
| 1297 | !---------------------------------------------------------------------- |
---|
| 1298 | ! ... Local variables |
---|
| 1299 | !---------------------------------------------------------------------- |
---|
| 1300 | REAL,DIMENSION(nx,ny,nz) :: p3 ! pressures at model layers (full levels) (Pa) |
---|
| 1301 | REAL,DIMENSION(nx,ny,nz) :: t3 ! temp. at model layer (full levels) (K) |
---|
| 1302 | INTEGER :: i,j,k,l ! lon, lat and vertical indices |
---|
| 1303 | REAL, DIMENSION(nx,ny) :: ptropd ! 2D field of tropopause pressures |
---|
| 1304 | INTEGER, DIMENSION(nx,ny) :: itropd ! 2D field of tropopause level indices |
---|
| 1305 | REAL, DIMENSION(nx,ny) :: itropdr ! 2D field of tropopause level indices |
---|
| 1306 | REAL, DIMENSION(PLON_GLO) :: itropr |
---|
| 1307 | REAL, DIMENSION(nz) :: p1,t1 ! 1D columns of pressures & temperatures |
---|
| 1308 | INTEGER :: jlat60S, jlat60N ! latitude indices for 60N and 60S |
---|
| 1309 | REAL, PARAMETER :: rinval=-999.00 ! invalid real number |
---|
| 1310 | INTEGER :: l300hPa ! index of layer containing 300 hPa level |
---|
| 1311 | INTEGER :: l120hPa ! index of layer containing 120 hPa level |
---|
| 1312 | REAL :: ptpm60,ptpmlat,ptpsum ! mean and integrated TP pressures |
---|
| 1313 | REAL :: rscale ! scaling factor to correct non-sense |
---|
| 1314 | REAL :: deltalat |
---|
| 1315 | REAL :: ptrop_glo(PLON_GLO) |
---|
| 1316 | INTEGER :: itrop_glo(PLON_GLO) |
---|
| 1317 | INTEGER :: nbbeg_loc,nbend_loc |
---|
| 1318 | |
---|
| 1319 | |
---|
| 1320 | |
---|
| 1321 | deltalat = 180./FLOAT(ny-1) |
---|
| 1322 | |
---|
| 1323 | !T and p to 2D grid |
---|
| 1324 | CALL grid1dto2d_glo(pmid,p3) |
---|
| 1325 | CALL grid1dTo2d_glo(temp,t3) |
---|
| 1326 | |
---|
| 1327 | !indices of lat=60N and 60S |
---|
| 1328 | jlat60N=FLOOR(30./deltalat)+1 |
---|
| 1329 | jlat60S=ny-jlat60N+1 |
---|
| 1330 | |
---|
| 1331 | !Calc index of model layer (typically) containing the 300 hPa level |
---|
| 1332 | !loop from surface to top |
---|
| 1333 | DO l=1,nz-1 |
---|
| 1334 | IF (0.5*(LOG(p3(nx/2,ny/2,l))+LOG(p3(nx/2,ny/2,l+1))).LT.LOG(30000.)) & |
---|
| 1335 | GOTO 100 |
---|
| 1336 | ENDDO |
---|
| 1337 | 100 CONTINUE |
---|
| 1338 | l300hPa=l |
---|
| 1339 | |
---|
| 1340 | !Calc index of model layer (typically) containing the 120 hPa level |
---|
| 1341 | !loop from surface to top |
---|
| 1342 | DO l=1,nz-1 |
---|
| 1343 | IF (0.5*(LOG(p3(nx/2,ny/2,l))+LOG(p3(nx/2,ny/2,l+1))).LT.LOG(12000.)) & |
---|
| 1344 | GOTO 101 |
---|
| 1345 | ENDDO |
---|
| 1346 | 101 CONTINUE |
---|
| 1347 | l120hPa=l |
---|
| 1348 | |
---|
| 1349 | !loop over southern hemisphere and get tropopause (loop starts at equator!!) |
---|
| 1350 | DO j=ny/2,ny |
---|
| 1351 | ptpsum=0. ! initialize sum of tropopause pressures |
---|
| 1352 | DO i=1,nx |
---|
| 1353 | !reverse order of levels because stattrop |
---|
| 1354 | !requires first index at top of atmosphere |
---|
| 1355 | DO l=1,nz |
---|
| 1356 | p1(l)=p3(i,j,nz+1-l) |
---|
| 1357 | t1(l)=t3(i,j,nz+1-l) |
---|
| 1358 | ENDDO |
---|
| 1359 | CALL stattrop(p1,t1,nz,ptropd(i,j),itropd(i,j)) |
---|
| 1360 | itropd(i,j)=nz-itropd(i,j)+1 ! reverse tropopause level index back |
---|
| 1361 | |
---|
| 1362 | !if no tropopause was found (happens only rarely) |
---|
| 1363 | IF (ptropd(i,j).EQ.rinval) THEN |
---|
| 1364 | IF (j.GT.jlat60S) THEN |
---|
| 1365 | ptropd(i,j)=30000. ! set to 300 hPa at high latitudes |
---|
| 1366 | itropd(i,j)=l300hPa ! index of layer containing 300 hPa level |
---|
| 1367 | ELSE |
---|
| 1368 | ptropd(i,j)=12000. ! set to 120 hPa at low latitudes |
---|
| 1369 | itropd(i,j)=l120hPa ! index of layer containing 120 hPa level |
---|
| 1370 | END IF |
---|
| 1371 | ENDIF |
---|
| 1372 | |
---|
| 1373 | !sum up tropopause levels to calculate mean |
---|
| 1374 | ptpsum=ptpsum+ptropd(i,j) |
---|
| 1375 | ENDDO |
---|
| 1376 | |
---|
| 1377 | !Correct nonsense at high latitudes (correction is only |
---|
| 1378 | !necessary during southern hemisphere winter from Apr until Sep) |
---|
| 1379 | IF (j.EQ.jlat60S) ptpm60=ptpsum/float(nx) |
---|
| 1380 | |
---|
| 1381 | IF (j.GT.jlat60S) THEN |
---|
| 1382 | ptpmlat=ptpsum/float(nx) |
---|
| 1383 | !correct ptropd if mean tropopause pressure at this latitude |
---|
| 1384 | !is lower than at 60S |
---|
| 1385 | IF (ptpmlat.LT.ptpm60) THEN |
---|
| 1386 | rscale=ptpm60/ptpmlat |
---|
| 1387 | DO i=1,nx |
---|
| 1388 | ptropd(i,j)=ptropd(i,j)*rscale |
---|
| 1389 | !find corresponding level indices (start loop at surface) |
---|
| 1390 | DO l=1,nz-1 |
---|
| 1391 | IF (0.5*(LOG(p3(i,j,l))+LOG(p3(i,j,l+1))) & |
---|
| 1392 | .LT.LOG(ptropd(i,j))) GOTO 102 |
---|
| 1393 | ENDDO |
---|
| 1394 | 102 itropd(i,j)=l |
---|
| 1395 | ENDDO |
---|
| 1396 | ENDIF |
---|
| 1397 | ENDIF |
---|
| 1398 | ENDDO |
---|
| 1399 | |
---|
| 1400 | !loop over northern hemisphere and get tropopause (loop starts at equator!!) |
---|
| 1401 | DO j=ny/2,1,-1 |
---|
| 1402 | ptpsum=0. ! initialize sum of tropopause pressures |
---|
| 1403 | DO i=1,nx |
---|
| 1404 | !reverse order of levels because stattrop |
---|
| 1405 | !requires first index at top of atmosphere |
---|
| 1406 | DO l=1,nz |
---|
| 1407 | p1(l)=p3(i,j,nz+1-l) |
---|
| 1408 | t1(l)=t3(i,j,nz+1-l) |
---|
| 1409 | ENDDO |
---|
| 1410 | CALL stattrop(p1,t1,nz,ptropd(i,j),itropd(i,j)) |
---|
| 1411 | itropd(i,j)=nz-itropd(i,j)+1 ! reverse tropopause level index back |
---|
| 1412 | |
---|
| 1413 | !if no tropopause was found (happens only rarely) |
---|
| 1414 | IF (ptropd(i,j).EQ.rinval) THEN |
---|
| 1415 | IF (j.LT.jlat60N) THEN |
---|
| 1416 | ptropd(i,j)=30000. ! set to 300 hPa for high latitudes |
---|
| 1417 | itropd(i,j)=l300hPa ! index of layer containing 300 hPa level |
---|
| 1418 | ELSE |
---|
| 1419 | ptropd(i,j)=12000. ! set to 120 hPa for low latitudes |
---|
| 1420 | itropd(i,j)=l120hPa ! index of layer containing 120 hPa level |
---|
| 1421 | END IF |
---|
| 1422 | ENDIF |
---|
| 1423 | |
---|
| 1424 | !sum up tropopause levels to calculate mean |
---|
| 1425 | ptpsum=ptpsum+ptropd(i,j) |
---|
| 1426 | ENDDO |
---|
| 1427 | |
---|
| 1428 | !Correct nonsense at high latitudes (correction is only |
---|
| 1429 | !necessary during nouthern hemisphere winter from Oct-Mar) |
---|
| 1430 | IF (j.EQ.jlat60N) ptpm60=ptpsum/float(nx) |
---|
| 1431 | |
---|
| 1432 | IF (j.LT.jlat60N) THEN |
---|
| 1433 | ptpmlat=ptpsum/float(nx) |
---|
| 1434 | !correct ptropd if mean tropopause pressure at this latitude |
---|
| 1435 | !is lower than at 60S |
---|
| 1436 | IF (ptpmlat.LT.ptpm60) THEN |
---|
| 1437 | rscale=ptpm60/ptpmlat |
---|
| 1438 | DO i=1,nx |
---|
| 1439 | ptropd(i,j)=ptropd(i,j)*rscale |
---|
| 1440 | !find corresponding level indices (start loop at surface) |
---|
| 1441 | DO l=1,nz-1 |
---|
| 1442 | IF (0.5*(LOG(p3(i,j,l))+LOG(p3(i,j,l+1))) & |
---|
| 1443 | .LT.LOG(ptropd(i,j))) GOTO 103 |
---|
| 1444 | ENDDO |
---|
| 1445 | 103 itropd(i,j)=l |
---|
| 1446 | ENDDO |
---|
| 1447 | ENDIF |
---|
| 1448 | ENDIF |
---|
| 1449 | ENDDO |
---|
| 1450 | |
---|
| 1451 | ! borne pour le passage GLO -> LOC |
---|
| 1452 | nbbeg_loc = nbp_para_begin(mpi_rank)+plon_omp_begin-1 |
---|
| 1453 | nbend_loc = nbp_para_begin(mpi_rank)+plon_omp_end-1 |
---|
| 1454 | !Back to physiq grid 1D |
---|
| 1455 | CALL grid2dto1d_glo(ptropd,ptrop_glo) |
---|
| 1456 | ptrop=ptrop_glo(nbbeg_loc:nbend_loc) |
---|
| 1457 | |
---|
| 1458 | |
---|
| 1459 | itropdr = FLOAT(itropd) |
---|
| 1460 | CALL grid2dto1d_glo(itropdr,itropr) |
---|
| 1461 | itrop_glo = INT(itropr) |
---|
| 1462 | |
---|
| 1463 | itrop=itrop_glo(nbbeg_loc:nbend_loc) |
---|
| 1464 | |
---|
| 1465 | END SUBROUTINE FDTROPOPAUSE |
---|
| 1466 | |
---|
| 1467 | SUBROUTINE PINTERP(ncol, & |
---|
| 1468 | fieldin, & |
---|
| 1469 | pin, & |
---|
| 1470 | nin, & |
---|
| 1471 | fieldout, & |
---|
| 1472 | pout, & |
---|
| 1473 | nout, & |
---|
| 1474 | index, & |
---|
| 1475 | entered) |
---|
| 1476 | !----------------------------------------------------------------------- |
---|
| 1477 | ! ... Interpolates on model pressure levels |
---|
| 1478 | ! Stacy Walters and Didier Hauglustaine, 1999. |
---|
| 1479 | !----------------------------------------------------------------------- |
---|
| 1480 | |
---|
| 1481 | IMPLICIT NONE |
---|
| 1482 | |
---|
| 1483 | !----------------------------------------------------------------------- |
---|
| 1484 | ! ... Dummy args |
---|
| 1485 | !----------------------------------------------------------------------- |
---|
| 1486 | LOGICAL, INTENT(in) :: entered |
---|
| 1487 | INTEGER, INTENT(in) :: nin |
---|
| 1488 | INTEGER, INTENT(in) :: nout |
---|
| 1489 | INTEGER, INTENT(in) :: ncol |
---|
| 1490 | INTEGER, INTENT(inout) :: INDEX(ncol,nout) |
---|
| 1491 | REAL, INTENT(in) :: pin(ncol,nin) |
---|
| 1492 | REAL, INTENT(in) :: fieldin(ncol,nin) |
---|
| 1493 | REAL, INTENT(in) :: pout(ncol,nout) |
---|
| 1494 | REAL, INTENT(out) :: fieldout(ncol,nout) |
---|
| 1495 | |
---|
| 1496 | !----------------------------------------------------------------------- |
---|
| 1497 | ! ... Local variables |
---|
| 1498 | !----------------------------------------------------------------------- |
---|
| 1499 | INTEGER :: i, j, k |
---|
| 1500 | REAL :: deltp(ncol,nout) |
---|
| 1501 | REAL :: pinratio(ncol,nin) |
---|
| 1502 | |
---|
| 1503 | IF ( .NOT. entered) THEN |
---|
| 1504 | |
---|
| 1505 | DO i = 1, ncol |
---|
| 1506 | DO k = 1, nout |
---|
| 1507 | |
---|
| 1508 | DO j = nin-1, 1, -1 |
---|
| 1509 | IF (pout(i,k) <= pin(i,j)) THEN |
---|
| 1510 | IF (pout(i,k) > pin(i,j+1)) THEN |
---|
| 1511 | index(i,k) = j |
---|
| 1512 | ENDIF |
---|
| 1513 | ENDIF |
---|
| 1514 | ENDDO |
---|
| 1515 | |
---|
| 1516 | ENDDO |
---|
| 1517 | ENDDO |
---|
| 1518 | |
---|
| 1519 | ENDIF |
---|
| 1520 | |
---|
| 1521 | DO i = 1, ncol |
---|
| 1522 | DO j = 1, nin-1 |
---|
| 1523 | pinratio(i,j) = 1./LOG(pin(i,j)/pin(i,j+1)) |
---|
| 1524 | END DO |
---|
| 1525 | END DO |
---|
| 1526 | |
---|
| 1527 | |
---|
| 1528 | DO i = 1, ncol |
---|
| 1529 | DO k = 1, nout |
---|
| 1530 | |
---|
| 1531 | |
---|
| 1532 | IF (pout(i,k) <= pin(i,nin)) THEN |
---|
| 1533 | fieldout(i,k) = fieldin(i,nin) |
---|
| 1534 | ELSE IF (pout(i,k) >= pin(i,1)) THEN |
---|
| 1535 | fieldout(i,k) = fieldin(i,1) |
---|
| 1536 | ELSE |
---|
| 1537 | |
---|
| 1538 | deltp(i,k) = LOG(pout(i,k)/pin(i,index(i,k)+1)) * pinratio(i,index(i,k)) |
---|
| 1539 | |
---|
| 1540 | fieldout(i,k) = fieldin(i,index(i,k)+1) & |
---|
| 1541 | + deltp(i,k) * (fieldin(i,index(i,k))-fieldin(i,index(i,k)+1)) |
---|
| 1542 | |
---|
| 1543 | ENDIF |
---|
| 1544 | |
---|
| 1545 | ENDDO |
---|
| 1546 | ENDDO |
---|
| 1547 | |
---|
| 1548 | END SUBROUTINE PINTERP |
---|
| 1549 | |
---|
| 1550 | SUBROUTINE stattrop(pfull, tfull, nlev, ptropd, itropd) |
---|
| 1551 | |
---|
| 1552 | !********************************************************************************** |
---|
| 1553 | ! |
---|
| 1554 | ! input : p, t real(nlev) |
---|
| 1555 | ! input : nlev real |
---|
| 1556 | ! output: ptropd real |
---|
| 1557 | ! output: itropd integer |
---|
| 1558 | ! |
---|
| 1559 | ! programmed by Dominik Brunner V1.0 Aug 2000 |
---|
| 1560 | ! |
---|
| 1561 | ! built upon routine stattrop by Peter van Velthoven, |
---|
| 1562 | ! KNMI, The Netherlands |
---|
| 1563 | ! and on the ECHAM4/CHEM routine xttphwmo by |
---|
| 1564 | ! Thomas Reichle, Christine Land, B. Steil and R. Hein, DLR |
---|
| 1565 | ! |
---|
| 1566 | ! purpose |
---|
| 1567 | ! ------- |
---|
| 1568 | ! stattrop computes the pressure (Pa) at the thermal (static) tropopause |
---|
| 1569 | ! (TP) for a 1D column (sounding) of pressures and temperatures following |
---|
| 1570 | ! the definition of the height of the tropopause as postulated by WMO (1957). |
---|
| 1571 | ! |
---|
| 1572 | ! ATTENTION: In the current formulation of the program the first |
---|
| 1573 | ! level (index 1) must be at the top of the atmosphere and the |
---|
| 1574 | ! last level (index nlev) at the surface |
---|
| 1575 | ! |
---|
| 1576 | ! interface |
---|
| 1577 | ! --------- |
---|
| 1578 | ! call stattrop(pfull, tfull, nlev, ptropd, itropd) |
---|
| 1579 | ! - Input |
---|
| 1580 | ! nlev : number of vertical levels |
---|
| 1581 | ! pfull: pressure in 1D column at nlev full levels (layers) |
---|
| 1582 | ! tfull: temperature in 1D column at nlev full levels |
---|
| 1583 | ! - Output |
---|
| 1584 | ! ptropd: height of the tropopause in Pa |
---|
| 1585 | ! itropd: index of layer containing the tropopause |
---|
| 1586 | ! |
---|
| 1587 | ! methods |
---|
| 1588 | ! ------- |
---|
| 1589 | ! - Lapse rate gamma = -dT/dz |
---|
| 1590 | ! Using the hydrostatic approximation |
---|
| 1591 | ! |
---|
| 1592 | ! dz = -dp/(g*rho) = -dp/p * R/g * T = -dlnp * R/g * T |
---|
| 1593 | ! |
---|
| 1594 | ! we get -dT/dz = dT/T * g/R *1/dlnp = dlnT/dlnp |
---|
| 1595 | ! |
---|
| 1596 | ! - Variables are assumed to vary linearly in log(pressure) |
---|
| 1597 | ! - The tropopause is the lowest level above 450 hPa fullfilling |
---|
| 1598 | ! the WMO criterium. If ptropd is < 85 hPa it is set to 85 hPa. |
---|
| 1599 | ! If no tropopause is found ptropd is set to -999. |
---|
| 1600 | ! |
---|
| 1601 | ! references |
---|
| 1602 | ! ---------- |
---|
| 1603 | ! |
---|
| 1604 | ! - WMO (1992): International meteorological vocabulary, Genf, 784pp.: |
---|
| 1605 | ! |
---|
| 1606 | ! 1. The first tropopause is defined as the lowest level at which |
---|
| 1607 | ! the lapse rate decreases to 2 deg C per kilometer or less, |
---|
| 1608 | ! provided also the average lapse rate between this level and |
---|
| 1609 | ! all higher levels within 2 kilometers does not exceed 2 deg C |
---|
| 1610 | ! |
---|
| 1611 | ! - Randel WJ, Wu F, Gaffen DJ, Interannual variability of the tropical |
---|
| 1612 | ! tropopause derived from radiosonde data and NCEP reanalyses, |
---|
| 1613 | ! JOURNAL OF GEOPHYSICAL RESEARCH, 105, 15509-15523, 2000. |
---|
| 1614 | ! |
---|
| 1615 | ! The following webpage clearifies the calculation of the tropopause |
---|
| 1616 | ! in the NCEP reanalysis: http://dss.ucar.edu/pub/reanalysis/FAQ.html |
---|
| 1617 | ! |
---|
| 1618 | ! For comparison NCEP reanalysis tropopause pressures can be obtained |
---|
| 1619 | ! from http://www.cdc.noaa.gov/cdc/reanalysis/reanalysis.shtml |
---|
| 1620 | ! |
---|
| 1621 | !********************************************************************************** |
---|
| 1622 | USE INCA_DIM |
---|
| 1623 | |
---|
| 1624 | IMPLICIT NONE |
---|
| 1625 | INTEGER nlev |
---|
| 1626 | REAL pfull(nlev), tfull(nlev) |
---|
| 1627 | REAL ptropd, p2km, lapseavg, lapsesum |
---|
| 1628 | INTEGER ilev, ilev1, itropd, count, l |
---|
| 1629 | |
---|
| 1630 | REAL lapse(PLEV), phalf(PLEV) |
---|
| 1631 | REAL grav, rgas, const, lapsec, rinval |
---|
| 1632 | PARAMETER (grav=9.80665, rgas=287.04, & |
---|
| 1633 | const=1000.*grav/rgas, & |
---|
| 1634 | lapsec=2.0, rinval=-999.0 ) |
---|
| 1635 | ! min. and max. allowed tropopause pressure |
---|
| 1636 | REAL pmin, pmax |
---|
| 1637 | PARAMETER (pmin=8500.,pmax=45000.) |
---|
| 1638 | ! deltaz is layer thickness used for second tropopause criterium |
---|
| 1639 | REAL deltaz |
---|
| 1640 | PARAMETER (deltaz=2.) |
---|
| 1641 | LOGICAL found |
---|
| 1642 | REAL p1, p2 , weight |
---|
| 1643 | |
---|
| 1644 | ptropd = rinval |
---|
| 1645 | itropd = 0 |
---|
| 1646 | found = .FALSE. |
---|
| 1647 | |
---|
| 1648 | ! Loop from model top to surface and calculate lapse rate gamma=-dT/dz (K/km) |
---|
| 1649 | DO ilev = 1, nlev-1 |
---|
| 1650 | ilev1 = ilev + 1 |
---|
| 1651 | lapse(ilev) = LOG(tfull(ilev))- LOG(tfull(ilev1)) |
---|
| 1652 | lapse(ilev) = lapse(ilev) / & |
---|
| 1653 | (LOG(pfull(ilev))- LOG(pfull(ilev1))) |
---|
| 1654 | lapse(ilev) = const * lapse(ilev) |
---|
| 1655 | phalf(ilev) = (pfull(ilev) + pfull(ilev1))*0.5 |
---|
| 1656 | END DO |
---|
| 1657 | ! Loop from surface to top to find (lowest) tropopause |
---|
| 1658 | DO ilev = nlev-2, 1, -1 |
---|
| 1659 | ! **** 1st test: -dT/dz less than 2 K/km and pressure LT pmax? **** |
---|
| 1660 | !Modified 07/2001 -- D. Hauglustaine |
---|
| 1661 | IF (lapse(ilev).LT.lapsec.AND.pfull(ilev).LT.pmax.AND.ilev.GE.5) THEN |
---|
| 1662 | IF (.NOT.found) THEN |
---|
| 1663 | ! Interpolate tropopause pressure linearly in log(p) |
---|
| 1664 | p1 = LOG(phalf(ilev)) |
---|
| 1665 | p2 = LOG(phalf(ilev+1)) |
---|
| 1666 | !Modified 09/2001 -- L. Jourdain |
---|
| 1667 | IF ( (lapse(ilev).NE.lapse(ilev+1)) .AND. (lapse(ilev+1).GE.lapsec) ) THEN |
---|
| 1668 | weight = (lapsec - lapse(ilev)) / & |
---|
| 1669 | (lapse(ilev+1)-lapse(ilev)) |
---|
| 1670 | ptropd = EXP(p1 + weight * (p2-p1)) ! tropopause pressure |
---|
| 1671 | ELSE |
---|
| 1672 | ptropd = phalf(ilev) |
---|
| 1673 | END IF |
---|
| 1674 | ! *** 2nd test: average -dT/dz in layer deltaz above TP must not be greater than 2K/km |
---|
| 1675 | p2km = EXP(LOG(ptropd)-deltaz*const/tfull(ilev)) ! press 2 km above TP |
---|
| 1676 | lapsesum = 0.0 ! init avg. lapse rate 2 km above TP |
---|
| 1677 | lapseavg=1.e9 |
---|
| 1678 | count = 0 |
---|
| 1679 | DO l=ilev,1,-1 |
---|
| 1680 | IF (phalf(l).LT.p2km) GOTO 100 |
---|
| 1681 | lapsesum=lapsesum+lapse(l) |
---|
| 1682 | count=count+1 |
---|
| 1683 | lapseavg=lapsesum/float(count) |
---|
| 1684 | END DO |
---|
| 1685 | 100 CONTINUE |
---|
| 1686 | found=lapseavg.LT.lapsec |
---|
| 1687 | ! Discard tropopause? |
---|
| 1688 | IF (.NOT.found) THEN |
---|
| 1689 | ptropd=rinval |
---|
| 1690 | ELSE |
---|
| 1691 | ptropd=MAX(ptropd,pmin) ! ptropd must be GE 85 hPa |
---|
| 1692 | itropd=ilev ! assign level index |
---|
| 1693 | END IF |
---|
| 1694 | END IF |
---|
| 1695 | END IF |
---|
| 1696 | END DO |
---|
| 1697 | RETURN |
---|
| 1698 | END SUBROUTINE stattrop |
---|
| 1699 | |
---|
| 1700 | |
---|
| 1701 | |
---|
| 1702 | subroutine CALC_PV(latitude_deg,paprs,pplay,ts,t,rot) |
---|
| 1703 | ! This subroutine consists in deriving Ertel's potential vorticity (PV) in |
---|
| 1704 | ! potential vorticity units (PVU), i.e. 1 PVU = 1E-6 K.m^2/kg/s. |
---|
| 1705 | ! It accounts for the vertical potential temperature gradient (=static stability) |
---|
| 1706 | ! and the (almost) horizontal wind vorticity. |
---|
| 1707 | ! The equation and clear explanations are available in Gettelman et al. (2011, Rev. Geophys), |
---|
| 1708 | ! a review paper on the extratropical UTLS. |
---|
| 1709 | ! Added by Yann Cohen, November 2019. ! Created by Yann Cohen |
---|
| 1710 | |
---|
| 1711 | USE CONST_MOD |
---|
| 1712 | USE XIOS_INCA |
---|
| 1713 | IMPLICIT NONE |
---|
| 1714 | |
---|
| 1715 | !------------------------------------------------------------------------------- |
---|
| 1716 | ! Arguments: |
---|
| 1717 | REAL, INTENT(IN) :: latitude_deg(PLON) ! latitude |
---|
| 1718 | REAL, INTENT(IN) :: paprs(PLON,PLEVP) !--- Cells-edges pressure |
---|
| 1719 | REAL, INTENT(IN) :: pplay(PLON,PLEV) !--- Cells-centers pressure |
---|
| 1720 | REAL, INTENT(IN) :: ts(PLON) !--- Surface temperature |
---|
| 1721 | REAL, INTENT(IN) :: t(PLON,PLEV) !--- Cells-centers temperature |
---|
| 1722 | REAL, INTENT(IN) :: rot(PLON,PLEV) !--- Cells-centers relative vorticity |
---|
| 1723 | !------------------------------------------------------------------------------- |
---|
| 1724 | ! Local variables: |
---|
| 1725 | include "YOMCST_I.h" |
---|
| 1726 | INTEGER :: i, k |
---|
| 1727 | REAL :: al,al2top |
---|
| 1728 | REAL,PARAMETER :: preff=101325 ! Unit: Pa. |
---|
| 1729 | REAL, DIMENSION(PLON,PLEVP):: tpot_edg |
---|
| 1730 | REAL, DIMENSION(PLON,PLEV) :: tpot_cen |
---|
| 1731 | REAL, DIMENSION(PLON,PLEV) :: PV |
---|
| 1732 | !------------------------------------------------------------------------------- |
---|
| 1733 | PV(:,:) = 0. |
---|
| 1734 | |
---|
| 1735 | !--- POTENTIAL TEMPERATURE AT CELLS CENTERS AND INTERFACES |
---|
| 1736 | DO i = 1,PLON |
---|
| 1737 | tpot_cen(i,1) = t(i,1)*(preff/pplay(i,1))**RKAPPA |
---|
| 1738 | tpot_edg(i,1) = ts(i) *(preff/paprs(i,1))**RKAPPA |
---|
| 1739 | |
---|
| 1740 | DO k=2,PLEV |
---|
| 1741 | |
---|
| 1742 | tpot_cen(i,k) = t(i,k)*(preff/pplay(i,k))**RKAPPA |
---|
| 1743 | al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k)) |
---|
| 1744 | tpot_edg(i,k) = (t(i,k-1)+al*(t(i,k)-t(i,k-1)))*(preff/paprs(i,k))**RKAPPA |
---|
| 1745 | |
---|
| 1746 | END DO |
---|
| 1747 | |
---|
| 1748 | ! Last ingredient for the PV field on the whole vertical grid: |
---|
| 1749 | ! linear extrapolation up to the summit, for cell-edges potential temperature. |
---|
| 1750 | ! The subsequent Theta vertical gradient is then used for the PV calculation at the top (full) level. |
---|
| 1751 | ! al2top=LOG(paprs(i,PLEV)/paprs(i,PLEVP))/LOG(paprs(i,PLEV)/pplay(i,PLEV)) |
---|
| 1752 | ! tpot_edg(i,PLEVP) = tpot_edg(i,PLEV) + al2top*(tpot_cen(i,PLEV)-tpot_edg(i,PLEV)) |
---|
| 1753 | |
---|
| 1754 | END DO |
---|
| 1755 | |
---|
| 1756 | |
---|
| 1757 | !--- ERTEL POTENTIAL VORTICITY AT CELLS CENTERS (except in top layer) |
---|
| 1758 | DO i = 1, PLON |
---|
| 1759 | DO k= 1, PLEV-1 |
---|
| 1760 | PV(i,k)=-1.E6*RG*(rot(i,k)+2.*ROMEGA*SIN(latitude_deg(i)*RPI/180.))& |
---|
| 1761 | * (tpot_edg(i,k+1)-tpot_edg(i,k)) / (paprs(i,k+1)-paprs(i,k)) |
---|
| 1762 | END DO |
---|
| 1763 | END DO |
---|
| 1764 | |
---|
| 1765 | |
---|
| 1766 | CALL xios_inca_change_context("inca") |
---|
| 1767 | CALL xios_inca_send_field("PV", PV) |
---|
| 1768 | CALL xios_inca_change_context("LMDZ") |
---|
| 1769 | |
---|
| 1770 | end subroutine CALC_PV |
---|