Ignore:
Timestamp:
Nov 1, 2016, 11:19:45 AM (8 years ago)
Author:
oboucher
Message:

Removing useless global variables
Interpolation only done once a month
Changing test orders in implicit coagulation routine as some compilers didn't like it

Location:
LMDZ5/trunk/libf/phylmd/StratAer
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/StratAer/cond_evap_tstep_mod.F90

    r2690 r2695  
    1818      USE aerophys
    1919      USE infotrac
    20       USE YOMCST
     20      USE YOMCST, ONLY : RPI
    2121
    2222      IMPLICIT NONE
    … …  
    139139      USE aerophys
    140140      USE infotrac
     141      USE YOMCST, ONLY : RPI
    141142
    142143      IMPLICIT NONE
    143 
    144       include "YOMCST.h"
    145144
    146145      ! input variables
    … …  
    175174        ff(IK,:)=0.0
    176175        DO k=1, nbtr_bin
    177           IF (k.LE.(nbtr_bin-1).AND.Vbin_wet(k).LE.Vnew.AND.Vnew.LT.Vbin_wet(k+1)) THEN
    178             ff(IK,k)= Vbin_wet(k)/Vnew*(Vbin_wet(k+1)-Vnew)/(Vbin_wet(k+1)-Vbin_wet(k))
     176          IF (k.LE.(nbtr_bin-1)) THEN
     177            IF (Vbin_wet(k).LE.Vnew.AND.Vnew.LT.Vbin_wet(k+1)) THEN
     178              ff(IK,k)= Vbin_wet(k)/Vnew*(Vbin_wet(k+1)-Vnew)/(Vbin_wet(k+1)-Vbin_wet(k))
     179            ENDIF
    179180          ENDIF
    180181          IF (k.EQ.1.AND.Vnew.LE.Vbin_wet(k)) THEN
    181182            ff(IK,k)= 1.
    182183          ENDIF
    183           IF (k.GT.1.AND.Vbin_wet(k-1).LT.Vnew.AND.Vnew.LT.Vbin_wet(k)) THEN
    184             ff(IK,k)= 1.-ff(IK,k-1)
     184          IF (k.GT.1) THEN
     185            IF (Vbin_wet(k-1).LT.Vnew.AND.Vnew.LT.Vbin_wet(k)) THEN
     186              ff(IK,k)= 1.-ff(IK,k-1)
     187            ENDIF
    185188          ENDIF
    186189          IF (k.EQ.nbtr_bin.AND.Vnew.GE.Vbin_wet(k)) THEN

  • LMDZ5/trunk/libf/phylmd/StratAer/interp_sulf_input.F90

    r2690 r2695  
    1 SUBROUTINE interp_sulf_input(debutphy,pdtphys,paprs,tr_seri,SO2_backgr_tend, &
    2                             & OCS_backgr_tend,SO2_lifetime,OCS_lifetime)
     1SUBROUTINE interp_sulf_input(debutphy,pdtphys,paprs,tr_seri)
    32
    43  USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_dimid, &
    … …  
    87  USE mod_grid_phy_lmdz
    98  USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
     9  USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
     10  USE phys_local_var_mod, ONLY : OCS_backgr_tend, SO2_backgr_tend
     11  USE phys_local_var_mod, ONLY : OCS_lifetime, SO2_lifetime
    1012  USE mod_phys_lmdz_para
    1113  USE dimphy
    … …  
    1315  USE infotrac
    1416  USE aerophys
     17  USE YOMCST
    1518
    1619  IMPLICIT NONE
    1720
    18   include "YOMCST.h"
    1921  include "dimensions.h"
    2022
    … …  
    2628
    2729! Variables locales
    28   INTEGER nmth
    2930  INTEGER n_lat   ! number of latitudes in the input data
    3031  INTEGER n_lon   ! number of longitudes in the input data
    3132  INTEGER, SAVE :: n_lev   ! number of levels in the input data
    3233  INTEGER n_mth   ! number of months in the input data
    33   REAL OCS_tmp
    34   REAL SO2_tmp
     34  REAL OCS_tmp, SO2_tmp
     35  INTEGER, SAVE :: mth_pre
    3536
    3637! Champs reconstitues
    3738  REAL paprs_glo(klon_glo,klev+1)
    38   REAL tr_seri_glo(klon_glo,klev,nbtr)
    3939
    4040  REAL, POINTER:: latitude(:)
    … …  
    5050! levels of input data
    5151
    52   REAL, ALLOCATABLE :: OCS_input(:, :, :, :)
    53   REAL, ALLOCATABLE :: SO2_input(:, :, :, :)
    54   REAL, ALLOCATABLE :: OCS_input_mth(:, :, :)
    55   REAL, ALLOCATABLE :: SO2_input_mth(:, :, :)
    56   REAL, SAVE, ALLOCATABLE :: OCS_input_tmp(:, :)
    57   REAL, SAVE, ALLOCATABLE :: SO2_input_tmp(:, :)
     52  REAL, ALLOCATABLE :: OCS_clim_in(:, :, :, :)
     53  REAL, ALLOCATABLE :: SO2_clim_in(:, :, :, :)
     54  REAL, ALLOCATABLE :: OCS_clim_mth(:, :, :)
     55  REAL, ALLOCATABLE :: SO2_clim_mth(:, :, :)
     56  REAL, ALLOCATABLE :: OCS_clim_tmp(:, :)
     57  REAL, ALLOCATABLE :: SO2_clim_tmp(:, :)
     58  REAL OCS_clim_glo(klon_glo,klev)
     59  REAL SO2_clim_glo(klon_glo,klev)
    5860  REAL, ALLOCATABLE :: OCS_lifetime_in(:, :, :, :)
    5961  REAL, ALLOCATABLE :: SO2_lifetime_in(:, :, :, :)
    6062  REAL, ALLOCATABLE :: OCS_lifetime_mth(:, :, :)
    6163  REAL, ALLOCATABLE :: SO2_lifetime_mth(:, :, :)
    62   REAL, SAVE, ALLOCATABLE :: OCS_lifetime_tmp(:, :)
    63   REAL, SAVE, ALLOCATABLE :: SO2_lifetime_tmp(:, :)
     64  REAL, ALLOCATABLE :: OCS_lifetime_tmp(:, :)
     65  REAL, ALLOCATABLE :: SO2_lifetime_tmp(:, :)
     66  REAL OCS_lifetime_glo(klon_glo,klev)
     67  REAL SO2_lifetime_glo(klon_glo,klev)
     68!
     69  REAL, ALLOCATABLE, SAVE :: OCS_clim(:,:)
     70  REAL, ALLOCATABLE, SAVE :: SO2_clim(:,:)
     71!$OMP THREADPRIVATE(OCS_clim,SO2_clim)
    6472!
    6573  INTEGER i, k, kk, ilon, ilev, j
    … …  
    6977  INTEGER ncid_in  ! IDs for input files
    7078  INTEGER varid, ncerr
    71 
    72 ! variable output
    73   REAL OCS_backgr_tend_glo(klon_glo,klev)
    74   REAL SO2_backgr_tend_glo(klon_glo,klev)
    75   REAL OCS_backgr_tend(klon,klev)
    76   REAL SO2_backgr_tend(klon,klev)
    77   REAL OCS_lifetime_glo(klon_glo,klev)
    78   REAL SO2_lifetime_glo(klon_glo,klev)
    79   REAL OCS_lifetime(klon,klev)
    80   REAL SO2_lifetime(klon,klev)
    8179   
    8280  INTEGER, PARAMETER :: lev_input=17
    … …  
    9694  REAL, PARAMETER :: min_SO2_lifetime=86400. !minimum SO2 lifetime [sec]
    9795
    98 !--preparation of fields
     96 IF (.NOT.ALLOCATED(OCS_clim)) ALLOCATE(OCS_clim(klon,klev))
     97 IF (.NOT.ALLOCATED(SO2_clim)) ALLOCATE(SO2_clim(klon,klev))
     98
     99  IF (debutphy.OR.mth_cur.NE.mth_pre) THEN
     100
     101!--preparation of global fields
    99102  CALL gather(paprs, paprs_glo)
    100   CALL gather(tr_seri, tr_seri_glo)
    101 
    102   IF (debutphy.AND.is_mpi_root) THEN
     103
     104  IF (is_mpi_root.AND.is_omp_root) THEN
    103105
    104106!--reading emission files
    … …  
    121123    n_mth = size(time)
    122124
    123     IF (.NOT.allocated(OCS_input))  allocate(OCS_input(n_lon, n_lat, n_lev, n_mth))
    124     IF (.NOT.allocated(SO2_input))  allocate(SO2_input(n_lon, n_lat, n_lev, n_mth))
    125     IF (.NOT.allocated(OCS_lifetime_in))  allocate(OCS_lifetime_in(n_lon, n_lat, n_lev, n_mth))
    126     IF (.NOT.allocated(SO2_lifetime_in))  allocate(SO2_lifetime_in(n_lon, n_lat, n_lev, n_mth))
     125    IF (.NOT.ALLOCATED(OCS_clim_in))     ALLOCATE(OCS_clim_in(n_lon, n_lat, n_lev, n_mth))
     126    IF (.NOT.ALLOCATED(SO2_clim_in))     ALLOCATE(SO2_clim_in(n_lon, n_lat, n_lev, n_mth))
     127    IF (.NOT.ALLOCATED(OCS_lifetime_in)) ALLOCATE(OCS_lifetime_in(n_lon, n_lat, n_lev, n_mth))
     128    IF (.NOT.ALLOCATED(SO2_lifetime_in)) ALLOCATE(SO2_lifetime_in(n_lon, n_lat, n_lev, n_mth))
    127129
    128130    CALL nf95_inq_varid(ncid_in, "OCS", varid)
    129     ncerr = nf90_get_var(ncid_in, varid, OCS_input)
     131    ncerr = nf90_get_var(ncid_in, varid, OCS_clim_in)
    130132    print *,'code erreur OCS=', ncerr, varid
    131133
    132134    CALL nf95_inq_varid(ncid_in, "SO2", varid)
    133     ncerr = nf90_get_var(ncid_in, varid, SO2_input)
     135    ncerr = nf90_get_var(ncid_in, varid, SO2_clim_in)
    134136    print *,'code erreur SO2=', ncerr, varid
    135137
    … …  
    144146    CALL nf95_close(ncid_in)
    145147
    146     IF (.NOT.allocated(OCS_input_mth)) allocate(OCS_input_mth(n_lon, n_lat, n_lev))
    147     IF (.NOT.allocated(SO2_input_mth)) allocate(SO2_input_mth(n_lon, n_lat, n_lev))
    148     IF (.NOT.allocated(OCS_input_tmp)) allocate(OCS_input_tmp(klon_glo, n_lev))
    149     IF (.NOT.allocated(SO2_input_tmp)) allocate(SO2_input_tmp(klon_glo, n_lev))
    150     IF (.NOT.allocated(OCS_lifetime_mth)) allocate(OCS_lifetime_mth(n_lon, n_lat, n_lev))
    151     IF (.NOT.allocated(SO2_lifetime_mth)) allocate(SO2_lifetime_mth(n_lon, n_lat, n_lev))
    152     IF (.NOT.allocated(OCS_lifetime_tmp)) allocate(OCS_lifetime_tmp(klon_glo, n_lev))
    153     IF (.NOT.allocated(SO2_lifetime_tmp)) allocate(SO2_lifetime_tmp(klon_glo, n_lev))
     148    IF (.NOT.ALLOCATED(OCS_clim_mth)) ALLOCATE(OCS_clim_mth(n_lon, n_lat, n_lev))
     149    IF (.NOT.ALLOCATED(SO2_clim_mth)) ALLOCATE(SO2_clim_mth(n_lon, n_lat, n_lev))
     150    IF (.NOT.ALLOCATED(OCS_clim_tmp)) ALLOCATE(OCS_clim_tmp(klon_glo, n_lev))
     151    IF (.NOT.ALLOCATED(SO2_clim_tmp)) ALLOCATE(SO2_clim_tmp(klon_glo, n_lev))
     152    IF (.NOT.ALLOCATED(OCS_lifetime_mth)) ALLOCATE(OCS_lifetime_mth(n_lon, n_lat, n_lev))
     153    IF (.NOT.ALLOCATED(SO2_lifetime_mth)) ALLOCATE(SO2_lifetime_mth(n_lon, n_lat, n_lev))
     154    IF (.NOT.ALLOCATED(OCS_lifetime_tmp)) ALLOCATE(OCS_lifetime_tmp(klon_glo, n_lev))
     155    IF (.NOT.ALLOCATED(SO2_lifetime_tmp)) ALLOCATE(SO2_lifetime_tmp(klon_glo, n_lev))
    154156
    155157!---select the correct month, undo multiplication with 1.e12 (precision reasons)
    156158!---correct latitudinal order and convert input from volume mixing ratio to mass mixing ratio
    157     nmth=mth_cur
    158159    DO j=1,n_lat
    159       SO2_input_mth(:,j,:) = 1.e-12*SO2_input(:,n_lat+1-j,:,nmth)*mSO2mol/mAIRmol
    160       OCS_input_mth(:,j,:) = 1.e-12*OCS_input(:,n_lat+1-j,:,nmth)*mOCSmol/mAIRmol
    161       SO2_lifetime_mth(:,j,:) = SO2_lifetime_in(:,n_lat+1-j,:,nmth)
    162       OCS_lifetime_mth(:,j,:) = OCS_lifetime_in(:,n_lat+1-j,:,nmth)
     160      SO2_clim_mth(:,j,:) = 1.e-12*SO2_clim_in(:,n_lat+1-j,:,mth_cur)*mSO2mol/mAIRmol
     161      OCS_clim_mth(:,j,:) = 1.e-12*OCS_clim_in(:,n_lat+1-j,:,mth_cur)*mOCSmol/mAIRmol
     162      SO2_lifetime_mth(:,j,:) = SO2_lifetime_in(:,n_lat+1-j,:,mth_cur)
     163      OCS_lifetime_mth(:,j,:) = OCS_lifetime_in(:,n_lat+1-j,:,mth_cur)
    163164    ENDDO
    164165
    165166!---reduce to a klon_glo grid but keep the levels
    166     CALL grid2dTo1d_glo(OCS_input_mth,OCS_input_tmp)
    167     CALL grid2dTo1d_glo(SO2_input_mth,SO2_input_tmp)
     167    CALL grid2dTo1d_glo(OCS_clim_mth,OCS_clim_tmp)
     168    CALL grid2dTo1d_glo(SO2_clim_mth,SO2_clim_tmp)
    168169    CALL grid2dTo1d_glo(OCS_lifetime_mth,OCS_lifetime_tmp)
    169170    CALL grid2dTo1d_glo(SO2_lifetime_mth,SO2_lifetime_tmp)
    … …  
    181182    ENDDO
    182183
    183   ENDIF !(debutphy.AND.is_mpi_root)
    184 
    185   IF (is_mpi_root) THEN
    186 
    187 !--set to background value everywhere in the very beginning, later only in the troposphere
    188       IF (debutphy.AND.MAXVAL(tr_seri_glo).LT.1.e-30) THEN
    189         p_bound=0.0
    190       ELSE
    191         p_bound=50000.
    192       ENDIF
    193 
    194 !--regridding tracer concentration on the vertical
    195       DO i=1, klon_glo
    196         DO k=1, klev
    197           !
    198           OCS_tmp=tr_seri_glo(i,k,id_OCS_strat)
    199           SO2_tmp=tr_seri_glo(i,k,id_SO2_strat)
    200           !--OCS and SO2 prescribed below p_bound
    201           IF (paprs_glo(i,k).GT.p_bound) THEN
    202             tr_seri_glo(i,k,id_OCS_strat)=0.0
    203             tr_seri_glo(i,k,id_SO2_strat)=0.0
    204             DO kk=1, n_lev
    205             tr_seri_glo(i,k,id_OCS_strat)=tr_seri_glo(i,k,id_OCS_strat)+ &
    206                MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
    207                *OCS_input_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
    208             tr_seri_glo(i,k,id_SO2_strat)=tr_seri_glo(i,k,id_SO2_strat)+ &
    209                MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
    210                *SO2_input_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
    211             ENDDO
    212           ENDIF
    213           OCS_backgr_tend_glo(i,k)=tr_seri_glo(i,k,id_OCS_strat)-OCS_tmp
    214           SO2_backgr_tend_glo(i,k)=tr_seri_glo(i,k,id_SO2_strat)-SO2_tmp
    215           !---regrid weighted lifetime
    216           OCS_lifetime_glo(i,k)=0.0
    217           SO2_lifetime_glo(i,k)=0.0
    218           DO kk=1, n_lev
    219           OCS_lifetime_glo(i,k)=OCS_lifetime_glo(i,k)+ &
    220                MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
    221                *OCS_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
    222           SO2_lifetime_glo(i,k)=SO2_lifetime_glo(i,k)+ &
    223                MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
    224                *SO2_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
    225           ENDDO
    226         ENDDO
     184  !---regrid weighted lifetime and climatologies
     185  DO i=1, klon
     186    DO k=1, klev
     187     OCS_lifetime_glo(i,k)=0.0
     188     SO2_lifetime_glo(i,k)=0.0
     189     OCS_clim_glo(i,k)=0.0
     190     SO2_clim_glo(i,k)=0.0
     191     DO kk=1, n_lev
     192      OCS_lifetime_glo(i,k)=OCS_lifetime_glo(i,k)+ &
     193           MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
     194           *OCS_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
     195      SO2_lifetime_glo(i,k)=SO2_lifetime_glo(i,k)+ &
     196           MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
     197           *SO2_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
     198      OCS_clim_glo(i,k)=OCS_clim_glo(i,k)+ &
     199           MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
     200           *OCS_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
     201      SO2_clim_glo(i,k)=SO2_clim_glo(i,k)+ &
     202           MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
     203           *SO2_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
    227204      ENDDO
    228 
    229   ENDIF !--is_mpi_root
    230 
    231   CALL scatter(tr_seri_glo, tr_seri)
    232   CALL scatter(OCS_backgr_tend_glo, OCS_backgr_tend)
    233   CALL scatter(SO2_backgr_tend_glo, SO2_backgr_tend)
     205    ENDDO
     206  ENDDO
     207
     208  ENDIF ! is_mpi_root
     209
     210!--keep memory of previous month
     211  mth_pre=mth_cur
     212
     213!--scatter global fields around
     214  CALL scatter(OCS_clim_glo, OCS_clim)
     215  CALL scatter(SO2_clim_glo, SO2_clim)
    234216  CALL scatter(OCS_lifetime_glo, OCS_lifetime)
    235217  CALL scatter(SO2_lifetime_glo, SO2_lifetime)
     218
     219  IF (is_mpi_root.AND.is_omp_root) THEN
     220!
     221    DEALLOCATE(OCS_clim_in,SO2_clim_in)
     222    DEALLOCATE(OCS_clim_mth,SO2_clim_mth)
     223    DEALLOCATE(OCS_clim_tmp,SO2_clim_tmp)
     224    DEALLOCATE(OCS_lifetime_in,SO2_lifetime_in)
     225    DEALLOCATE(OCS_lifetime_mth,SO2_lifetime_mth)
     226    DEALLOCATE(OCS_lifetime_tmp,SO2_lifetime_tmp)
     227!
     228  ENDIF !--is_mpi_root and is_omp_root
     229
     230  ENDIF ! debutphy.OR.new month
     231
     232!--set to background value everywhere in the very beginning, later only in the troposphere
     233!--a little dangerous as the MAXVAL is not computed on the global field
     234  IF (debutphy.AND.MAXVAL(tr_seri).LT.1.e-30) THEN
     235    p_bound=0.0
     236  ELSE
     237    p_bound=50000.
     238  ENDIF
     239
     240!--regridding tracer concentration on the vertical
     241  DO i=1, klon
     242    DO k=1, klev
     243      !
     244      OCS_tmp=tr_seri(i,k,id_OCS_strat)
     245      SO2_tmp=tr_seri(i,k,id_SO2_strat)
     246      !--OCS and SO2 prescribed below p_bound
     247      IF (paprs(i,k).GT.p_bound) THEN
     248        tr_seri(i,k,id_OCS_strat)=OCS_clim(i,k)
     249        tr_seri(i,k,id_SO2_strat)=SO2_clim(i,k)
     250      ENDIF
     251      OCS_backgr_tend(i,k)=tr_seri(i,k,id_OCS_strat)-OCS_tmp
     252      SO2_backgr_tend(i,k)=tr_seri(i,k,id_SO2_strat)-SO2_tmp
     253    ENDDO
     254  ENDDO
    236255
    237256  !convert SO2_backgr_tend from kg(SO2)/kgA to kg(S)/m2/layer/s for saving as diagnostic

  • LMDZ5/trunk/libf/phylmd/StratAer/micphy_tstep.F90

    r2690 r2695  
    88  USE cond_evap_tstep_mod
    99  USE sulfate_aer_mod, ONLY : STRAACT
    10   USE YOMCST
     10  USE YOMCST, ONLY : RPI, RD, RG
    1111
    1212  IMPLICIT NONE
    … …  
    7474    PDT=pdtphys
    7575    count_tstep=0
    76     DO while(PDT>0.0)
     76    DO WHILE (PDT>0.0)
    7777      count_tstep=count_tstep+1
    78       if( count_tstep .GT. nbtstep )  exit
     78      IF (count_tstep .GT. nbtstep)  EXIT
    7979      ! convert tr_seri(GASH2SO4) (in kg/kgA) to H2SO4 number density (in molecules/cm3)
    8080      rhoa=tr_seri(ilon,ilev,id_H2SO4_strat) &
    … …  
    145145  ENDDO
    146146
    147   IF (MINVAL(tr_seri(:,:,:)).LT.0.0) THEN
     147  IF (MINVAL(tr_seri).LT.0.0) THEN
    148148    DO ilon=1, klon
    149149    DO ilev=1, klev   

  • LMDZ5/trunk/libf/phylmd/StratAer/nucleation_tstep_mod.F90

    r2690 r2695  
    77  USE aerophys
    88  USE infotrac
    9   USE YOMCST
     9  USE YOMCST, ONLY : RPI, RD
    1010
    1111  IMPLICIT NONE
    … …  
    5252  USE aerophys
    5353  USE infotrac
    54   USE YOMCST
    5554
    5655  IMPLICIT NONE
    … …  
    7877  DO k=1, nbtr_bin
    7978  ! CK 20160531: bug fix for first bin
    80     IF (k.LE.(nbtr_bin-1).AND.Vbin(k).LE.Vnew.AND.Vnew.LT.Vbin(k+1)) THEN
    81       ff(k)= Vbin(k)/Vnew*(Vbin(k+1)-Vnew)/(Vbin(k+1)-Vbin(k))
     79    IF (k.LE.(nbtr_bin-1)) THEN
     80      IF (Vbin(k).LE.Vnew.AND.Vnew.LT.Vbin(k+1)) THEN
     81        ff(k)= Vbin(k)/Vnew*(Vbin(k+1)-Vnew)/(Vbin(k+1)-Vbin(k))
     82      ENDIF
    8283    ENDIF
    8384    IF (k.EQ.1.AND.Vnew.LE.Vbin(k)) THEN
    8485      ff(k)= 1.
    8586    ENDIF
    86     IF (k.GT.1.AND.Vbin(k-1).LT.Vnew.AND.Vnew.LT.Vbin(k)) THEN
    87       ff(k)= 1.-ff(k-1)
     87    IF (k.GT.1) THEN
     88      IF (Vbin(k-1).LT.Vnew.AND.Vnew.LT.Vbin(k)) THEN
     89        ff(k)= 1.-ff(k-1)
     90      ENDIF
    8891    ENDIF
    8992    IF (k.EQ.nbtr_bin.AND.Vnew.GE.Vbin(k)) THEN

  • LMDZ5/trunk/libf/phylmd/StratAer/ocs_to_so2.F90

    r2690 r2695  
    1 SUBROUTINE OCS_TO_SO2(pdtphys,tr_seri,t_seri,pplay,paprs,sh,OCS_lifetime,is_strato,ocs_convert_sub)
     1SUBROUTINE ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,sh,is_strato)
    22
    33  USE dimphy, ONLY : klon,klev
    … …  
    55  USE infotrac
    66  USE YOMCST, ONLY : RG
     7  USE phys_local_var_mod, ONLY : OCS_lifetime, ocs_convert
    78
    89  IMPLICIT NONE
    … …  
    1617  REAL,DIMENSION(klon,klev+1),INTENT(IN)        :: paprs   ! pression pour chaque inter-couche (en Pa)
    1718  REAL,DIMENSION(klon,klev),INTENT(IN)          :: sh      ! humidite specifique   
    18   REAL,DIMENSION(klon,klev),INTENT(IN)          :: OCS_lifetime
    1919  LOGICAL,DIMENSION(klon,klev),INTENT(IN)       :: is_strato
    2020
    2121  ! local variables
    2222  INTEGER                                       :: i,j,k,nb,ilon,ilev
    23   REAL,DIMENSION(klon,klev)                     :: ocs_convert_sub ! OCS converted to SO2 [kg(S)/m2/layer/s]
    2423
    2524!--convert OCS to SO2
     25  ocs_convert(:,:)=0.0
    2626  DO ilon=1, klon
    2727  DO ilev=1, klev
    … …  
    2929  IF (is_strato(ilon,ilev)) THEN
    3030    IF (OCS_lifetime(ilon,ilev).GT.0.0) THEN
    31       ocs_convert_sub(ilon,ilev)=tr_seri(ilon,ilev,id_OCS_strat)*(1-exp(-pdtphys/OCS_lifetime(ilon,ilev)))
    32     ELSE
    33       ocs_convert_sub(ilon,ilev)=0.0
     31      ocs_convert(ilon,ilev)=tr_seri(ilon,ilev,id_OCS_strat)*(1.0-exp(-pdtphys/OCS_lifetime(ilon,ilev)))
    3432    ENDIF
    35     tr_seri(ilon,ilev,id_OCS_strat)=tr_seri(ilon,ilev,id_OCS_strat) - ocs_convert_sub(ilon,ilev)
    36     tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) + mSO2mol/mOCSmol*ocs_convert_sub(ilon,ilev)
     33    tr_seri(ilon,ilev,id_OCS_strat)=tr_seri(ilon,ilev,id_OCS_strat) - ocs_convert(ilon,ilev)
     34    tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) + mSO2mol/mOCSmol*ocs_convert(ilon,ilev)
    3735    !convert ocs_convert from kg(OCS)/kgA to kg(S)/m2/layer/s for saving as diagnostic
    38     ocs_convert_sub(ilon,ilev)=ocs_convert_sub(ilon,ilev)*mSatom/mOCSmol* &
    39                             & (paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
     36    ocs_convert(ilon,ilev)=ocs_convert(ilon,ilev)*mSatom/mOCSmol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
    4037  ENDIF
    4138  ENDDO
    4239  ENDDO
    4340
    44 END SUBROUTINE OCS_TO_SO2
     41END SUBROUTINE ocs_to_so2

  • LMDZ5/trunk/libf/phylmd/StratAer/so2_to_h2so4.F90

    r2690 r2695  
    1 SUBROUTINE SO2_TO_H2SO4(pdtphys,tr_seri,t_seri,pplay,paprs,sh,SO2_lifetime,is_strato,sulf_convert_sub)
     1SUBROUTINE SO2_TO_H2SO4(pdtphys,tr_seri,t_seri,pplay,paprs,sh,is_strato)
    22
    33  USE dimphy, ONLY : klon,klev
    … …  
    55  USE infotrac
    66  USE YOMCST, ONLY : RG
     7  USE phys_local_var_mod, ONLY : SO2_lifetime, sulf_convert
    78
    89  IMPLICIT NONE
    … …  
    1718  REAL,DIMENSION(klon,klev+1),INTENT(IN)        :: paprs   ! pression pour chaque inter-couche (en Pa)
    1819  REAL,DIMENSION(klon,klev),INTENT(IN)          :: sh      ! humidite specifique   
    19   REAL,DIMENSION(klon,klev),INTENT(IN)          :: SO2_lifetime
    20   LOGICAL,DIMENSION(klon,klev),INTENT(IN)  :: is_strato
     20  LOGICAL,DIMENSION(klon,klev),INTENT(IN)       :: is_strato ! stratospheric flag
    2121
    2222  ! local variables in coagulation routine
    2323  INTEGER                                       :: i,j,k,nb,ilon,ilev
    24   REAL,DIMENSION(klon,klev)                     :: sulf_convert_sub ! SO2 converted to H2SO4 [kg(S)/m2/layer/s]
    2524
    2625!--convert SO2 to H2SO4
     26  sulf_convert(:,:)=0.0
    2727  DO ilon=1, klon
    2828  DO ilev=1, klev
    … …  
    3030  IF (is_strato(ilon,ilev)) THEN
    3131    IF (SO2_lifetime(ilon,ilev).GT.0.0) THEN
    32       sulf_convert_sub(ilon,ilev)=tr_seri(ilon,ilev,id_SO2_strat)*(1-exp(-pdtphys/SO2_lifetime(ilon,ilev)))
    33     ELSE
    34       sulf_convert_sub(ilon,ilev)=0.0
     32      sulf_convert(ilon,ilev)=tr_seri(ilon,ilev,id_SO2_strat)*(1.0-exp(-pdtphys/SO2_lifetime(ilon,ilev)))
    3533    ENDIF
    36     tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) - sulf_convert_sub(ilon,ilev)
    37     tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) + mH2SO4mol/mSO2mol*sulf_convert_sub(ilon,ilev)
     34    tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) - sulf_convert(ilon,ilev)
     35    tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) + mH2SO4mol/mSO2mol*sulf_convert(ilon,ilev)
    3836    !convert sulf_convert from kg(SO2)/kgA to kg(S)/m2/layer/s for saving as diagnostic
    39     sulf_convert_sub(ilon,ilev)=sulf_convert_sub(ilon,ilev)*mSatom/mSO2mol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
     37    sulf_convert(ilon,ilev)=sulf_convert(ilon,ilev)*mSatom/mSO2mol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
    4038  ENDIF
    4139  ENDDO

  • LMDZ5/trunk/libf/phylmd/StratAer/traccoag_mod.F90

    r2691 r2695  
    7676    REAL,PARAMETER    :: xlon_sai=120.35        ! longitude of SAI in degree
    7777
    78 !--be careful - this needs to be changed with resolution - here for 96x95
     78!--be careful - this needs to be changed with resolution - here for 96x95 - seems to work for 48x36 as well
    7979    REAL,PARAMETER    :: dlat=0.9474   ! d latitude in degree
    8080    REAL,PARAMETER    :: dlon=1.875    ! d longitude in degree
    … …  
    9898    REAL                                   :: zdz                 ! thickness of atm. model layer in m
    9999    REAL,DIMENSION(klon,klev)              :: dens_aer            ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
    100     REAL,DIMENSION(klon,klev)              :: sulf_convert_sub    ! SO2 converted to H2SO4 [kg(S)/m2/layer/s]
    101     REAL,DIMENSION(klon,klev)              :: sulf_nucl_sub       ! H2SO4 converted from gas to particles [kg(S)/m2/layer/s]
    102     REAL,DIMENSION(klon,klev)              :: ocs_convert_sub     ! OCS converted to SO2 [kg(S)/m2/layer/s]
    103     REAL,DIMENSION(klon,klev)              :: SO2_backgr_tend_sub ! SO2 from background [kg(S)/m2/layer/s]
    104     REAL,DIMENSION(klon,klev)              :: OCS_backgr_tend_sub ! OCS from background [kg(S)/m2/layer/s]
    105100
    106101    IF (is_mpi_root) THEN
    … …  
    134129!--initialising logical is_strato from stratomask
    135130    is_strato(:,:)=.FALSE.
    136     WHERE (stratomask(:,:).GT.0.5) is_strato(:,:)=.TRUE.
     131    WHERE (stratomask.GT.0.5) is_strato=.TRUE.
    137132
    138133! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
    … …  
    267262
    268263!--read background concentrations of OCS and SO2 and lifetimes from input file
    269     SO2_backgr_tend_sub(:,:)=0.0
    270     OCS_backgr_tend_sub(:,:)=0.0
    271     CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri,SO2_backgr_tend_sub, &
    272                           & OCS_backgr_tend_sub,SO2_lifetime,OCS_lifetime)
    273     SO2_backgr_tend(:,:)=SO2_backgr_tend_sub(:,:)
    274     OCS_backgr_tend(:,:)=OCS_backgr_tend_sub(:,:)
     264!--update the variables defined in phys_local_var_mod
     265    CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri)
    275266
    276267!--convert OCS to SO2 in the stratosphere
    277     ocs_convert_sub(:,:)=0.0
    278     CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,sh,OCS_lifetime,is_strato,ocs_convert_sub)
    279     ocs_convert(:,:)=ocs_convert_sub(:,:)
     268    CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,sh,is_strato)
    280269
    281270!--convert SO2 to H2SO4
    282     sulf_convert_sub(:,:)=0.0
    283     CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,sh,SO2_lifetime,is_strato,sulf_convert_sub)
    284     sulf_convert(:,:)=sulf_convert_sub(:,:)
     271    CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,sh,is_strato)
    285272
    286273!--common routine for nucleation and condensation/evaporation with adaptive timestep
    … …  
    307294    ENDDO
    308295
     296!    CALL minmaxsimple(tr_seri(:,:,id_SO2_strat),0.0,0.0,'fin SO2')
     297!    DO it=1, nbtr_bin
     298!      CALL minmaxsimple(tr_seri(:,:,nbtr_sulgas+it),0.0,0.0,'fin bin ')
     299!    ENDDO
     300
    309301  END SUBROUTINE traccoag
    310302
Note: See TracChangeset for help on using the changeset viewer.