Changeset 2788


Ignore:
Timestamp:
Feb 2, 2017, 7:01:50 PM (8 years ago)
Author:
dcugnet
Message:

Changes in ce0l about the way ozone forcing files are generated:

1) 3D raw input files "climoz.nc" are now handled.
2) Default behaviour is now to let the gcm interpolate in time online.

This helps to avoid huge forcing files (in particular for 3D fields).
In this case, the output files "climoz_LMDZ.nc" all have 14 records:

  • records 2-13 are obtained with records 1-12 of "climoz.nc".
  • records 1 and 14 are obtained respectively with:
    • record 12 of "climoz_m.nc" if available, of "climoz.nc" otherwise.
    • record 1 of "climoz_p.nc" if available, of "climoz.nc" otherwise.

3) If ok_daily_climoz key is TRUE, the time interpolation (one record

a day) is forced, using the 14 records described below.
This now depends on the calendar (it was on a 360 days basis only).

Changes in the gcm about the way zone forcing files are read/interpolated:

1) 3D horizontally interpolated "climoz_LMDZ.nc" files are now handled.
2) Daily files (already interpolated in time) are still handled, but their

number of records must match the expected number of days, that depends
on the calendar (records step is no longer 1/360 year).

3) 14 records monthly files are now handled (and prefered). This reduces

the I/O to a minimum and the aditional computational cost is low (simple
online linear time interpolation).

4) If adjust_tropopause key is TRUE, the input fields are stretched using

following method:

  • LMDZ dynamical tropopause is detected: Ptrop_lmdz = MAX ( P(Potential Vorticity==2PVU), P(theta==380K) )
  • file chemical tropopause is detected: Ptrop_file = P( tro3 == o3t ), where:

o3t = 91. + 28. * SIN(PI*(month-2)/6) (ppbV)

This formula comes from Thouret & al., ACP 6, 1033-1051, 2006.
The second term of the expression is multiplied by TANH(lat_deg/20.)
to account for latitude dependency.

  • File profile is streched in a +/- 5kms zone around the mean tropopause to ensure resulting tropopause matches the one of LMDZ. See procedure regr_pr_time_av for more details.
Location:
LMDZ5/trunk/libf
Files:
5 added
12 edited
5 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90

    r2738 r2788  
    8282  USE fonte_neige_mod
    8383  USE pbl_surface_mod
    84   USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz
     84  USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
    8585  USE indice_sol_mod
    8686  USE conf_phys_m, ONLY: conf_phys
    … …  
    170170  radsol(:) = 0.0                               !--- Net radiation at ground
    171171  rugmer(:) = 0.001                             !--- Ocean rugosity
    172   IF(read_climoz>=1) &                          !--- Ozone climatology
    173     CALL regr_lat_time_climoz(read_climoz)
     172  !--- Ozone (zonal or 3D) interpolation in space and time (if 2nd arg is TRUE)
     173  IF(read_climoz>=1) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
    174174
    175175! Sub-surfaces initialization.

  • LMDZ5/trunk/libf/misc/slopes_m.F90

    r2440 r2788  
    1 module slopes_m
     1MODULE slopes_m
    22
    33  ! Author: Lionel GUEZ
    4 
    5   implicit none
    6 
    7   interface slopes
    8      ! This generic function computes second order slopes with Van
    9      ! Leer slope-limiting, given cell averages. Reference: Dukowicz,
    10      ! 1987, SIAM Journal on Scientific and Statistical Computing, 8,
    11      ! 305.
    12 
    13      ! The only difference between the specific functions is the rank
    14      ! of the first argument and the equal rank of the result.
    15 
    16      ! real, intent(in), rank >= 1:: f ! (n, ...) cell averages, n must be >= 1
    17      ! real, intent(in):: x(:) ! (n + 1) cell edges
    18      ! real slopes, same shape as f ! (n, ...)
    19 
    20      module procedure slopes1, slopes2, slopes3, slopes4
    21   end interface
    22 
    23   private
    24   public slopes
    25 
    26 contains
    27 
    28   pure function slopes1(f, x)
    29 
    30     real, intent(in):: f(:)
    31     real, intent(in):: x(:)
    32     real slopes1(size(f))
    33 
    34     ! Local:
    35     integer n, i
    36     real xc(size(f)) ! (n) cell centers
    37     real h
    38 
    39     !------------------------------------------------------
    40 
    41     n = size(f)
    42     forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2.
    43     slopes1(1) = 0.
    44     slopes1(n) = 0.
    45 
    46     do i = 2, n - 1
    47        if (f(i) >= max(f(i - 1), f(i + 1)) .or. f(i) &
    48             <= min(f(i - 1), f(i + 1))) then
    49           ! Local extremum
    50           slopes1(i) = 0.
    51        else
    52           ! (f(i - 1), f(i), f(i + 1)) strictly monotonous
    53 
    54           ! Second order slope:
    55           slopes1(i) = (f(i + 1) - f(i - 1)) / (xc(i + 1) - xc(i - 1))
    56 
    57           ! Slope limitation:
    58           h = abs(x(i + 1) - xc(i))
    59           slopes1(i) = sign(min(abs(slopes1(i)), abs(f(i + 1) - f(i)) / h, &
    60                abs(f(i) - f(i - 1)) / h), slopes1(i))
    61        end if
    62     end do
    63 
    64   end function slopes1
    65 
    66   !*************************************************************
    67 
    68   pure function slopes2(f, x)
    69 
    70     real, intent(in):: f(:, :)
    71     real, intent(in):: x(:)
    72     real slopes2(size(f, 1), size(f, 2))
    73 
    74     ! Local:
    75     integer n, i, j
    76     real xc(size(f, 1)) ! (n) cell centers
    77     real h(2:size(f, 1) - 1), delta_xc(2:size(f, 1) - 1) ! (2:n - 1)
    78 
    79     !------------------------------------------------------
    80 
    81     n = size(f, 1)
    82     forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2.
    83 
    84     forall (i = 2:n - 1)
    85        h(i) = abs(x(i + 1) - xc(i))
    86        delta_xc(i) = xc(i + 1) - xc(i - 1)
    87     end forall
    88 
    89     do j = 1, size(f, 2)
    90        slopes2(1, j) = 0.
    91 
    92        do i = 2, n - 1
    93           if (f(i, j) >= max(f(i - 1, j), f(i + 1, j)) .or. &
    94                f(i, j) <= min(f(i - 1, j), f(i + 1, j))) then
    95              ! Local extremum
    96              slopes2(i, j) = 0.
    97           else
    98              ! (f(i - 1, j), f(i, j), f(i + 1, j))
    99              ! strictly monotonous
    100 
    101              ! Second order slope:
    102              slopes2(i, j) = (f(i + 1, j) - f(i - 1, j)) / delta_xc(i)
    103 
    104              ! Slope limitation:
    105              slopes2(i, j) = sign(min(abs(slopes2(i, j)), &
    106                   abs(f(i + 1, j) - f(i, j)) / h(i), &
    107                   abs(f(i, j) - f(i - 1, j)) / h(i)), slopes2(i, j))
    108           end if
    109        end do
    110 
    111        slopes2(n, j) = 0.
    112     end do
    113 
    114   end function slopes2
    115 
    116   !*************************************************************
    117 
    118   pure function slopes3(f, x)
    119 
    120     real, intent(in):: f(:, :, :)
    121     real, intent(in):: x(:)
    122     real slopes3(size(f, 1), size(f, 2), size(f, 3))
    123 
    124     ! Local:
    125     integer n, i, j, k
    126     real xc(size(f, 1)) ! (n) cell centers
    127     real h(2:size(f, 1) - 1), delta_xc(2:size(f, 1) - 1) ! (2:n - 1)
    128 
    129     !------------------------------------------------------
    130 
    131     n = size(f, 1)
    132     forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2.
    133 
    134     forall (i = 2:n - 1)
    135        h(i) = abs(x(i + 1) - xc(i))
    136        delta_xc(i) = xc(i + 1) - xc(i - 1)
    137     end forall
    138 
    139     do k = 1, size(f, 3)
    140        do j = 1, size(f, 2)
    141           slopes3(1, j, k) = 0.
    142 
    143           do i = 2, n - 1
    144              if (f(i, j, k) >= max(f(i - 1, j, k), f(i + 1, j, k)) .or. &
    145                   f(i, j, k) <= min(f(i - 1, j, k), f(i + 1, j, k))) then
    146                 ! Local extremum
    147                 slopes3(i, j, k) = 0.
    148              else
    149                 ! (f(i - 1, j, k), f(i, j, k), f(i + 1, j, k))
    150                 ! strictly monotonous
    151 
    152                 ! Second order slope:
    153                 slopes3(i, j, k) = (f(i + 1, j, k) - f(i - 1, j, k)) &
    154                      / delta_xc(i)
    155 
    156                 ! Slope limitation:
    157                 slopes3(i, j, k) = sign(min(abs(slopes3(i, j, k)), &
    158                      abs(f(i + 1, j, k) - f(i, j, k)) / h(i), &
    159                      abs(f(i, j, k) - f(i - 1, j, k)) / h(i)), slopes3(i, j, k))
    160              end if
    161           end do
    162 
    163           slopes3(n, j, k) = 0.
    164        end do
    165     end do
    166 
    167   end function slopes3
    168 
    169   !*************************************************************
    170 
    171   pure function slopes4(f, x)
    172 
    173     real, intent(in):: f(:, :, :, :)
    174     real, intent(in):: x(:)
    175     real slopes4(size(f, 1), size(f, 2), size(f, 3), size(f, 4))
    176 
    177     ! Local:
    178     integer n, i, j, k, l
    179     real xc(size(f, 1)) ! (n) cell centers
    180     real h(2:size(f, 1) - 1), delta_xc(2:size(f, 1) - 1) ! (2:n - 1)
    181 
    182     !------------------------------------------------------
    183 
    184     n = size(f, 1)
    185     forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2.
    186 
    187     forall (i = 2:n - 1)
    188        h(i) = abs(x(i + 1) - xc(i))
    189        delta_xc(i) = xc(i + 1) - xc(i - 1)
    190     end forall
    191 
    192     do l = 1, size(f, 4)
    193        do k = 1, size(f, 3)
    194           do j = 1, size(f, 2)
    195              slopes4(1, j, k, l) = 0.
    196 
    197              do i = 2, n - 1
    198                 if (f(i, j, k, l) >= max(f(i - 1, j, k, l), f(i + 1, j, k, l)) &
    199                      .or. f(i, j, k, l) &
    200                      <= min(f(i - 1, j, k, l), f(i + 1, j, k, l))) then
    201                    ! Local extremum
    202                    slopes4(i, j, k, l) = 0.
    203                 else
    204                    ! (f(i - 1, j, k, l), f(i, j, k, l), f(i + 1, j, k, l))
    205                    ! strictly monotonous
    206 
    207                    ! Second order slope:
    208                    slopes4(i, j, k, l) = (f(i + 1, j, k, l) &
    209                         - f(i - 1, j, k, l)) / delta_xc(i)
    210 
    211                    ! Slope limitation:
    212                    slopes4(i, j, k, l) = sign(min(abs(slopes4(i, j, k, l)), &
    213                         abs(f(i + 1, j, k, l) - f(i, j, k, l)) / h(i), &
    214                         abs(f(i, j, k, l) - f(i - 1, j, k, l)) / h(i)), &
    215                         slopes4(i, j, k, l))
    216                 end if
    217              end do
    218 
    219              slopes4(n, j, k, l) = 0.
    220           end do
    221        end do
    222     end do
    223 
    224   end function slopes4
    225 
    226 end module slopes_m
     4  ! Extension / factorisation: David CUGNET
     5
     6  IMPLICIT NONE
     7
     8  ! Those generic function computes second order slopes with Van
     9  ! Leer slope-limiting, given cell averages. Reference: Dukowicz,
     10  ! 1987, SIAM Journal on Scientific and Statistical Computing, 8,
     11  ! 305.
     12
     13  ! The only difference between the specific functions is the rank
     14  ! of the first argument and the equal rank of the result.
     15
     16  ! slope(ix,...) acts on ix th dimension.
     17
     18  ! real, intent(in), rank >= 1:: f ! (n, ...) cell averages, n must be >= 1
     19  ! real, intent(in):: x(:) ! (n + 1) cell edges
     20  ! real slopes, same shape as f ! (n, ...)
     21  INTERFACE slopes
     22     MODULE procedure slopes1, slopes2, slopes3, slopes4, slopes5
     23  END INTERFACE
     24
     25  PRIVATE
     26  PUBLIC :: slopes
     27
     28CONTAINS
     29
     30!-------------------------------------------------------------------------------
     31!
     32PURE FUNCTION slopes1(ix, f, x)
     33!
     34!-------------------------------------------------------------------------------
     35! Arguments:
     36  INTEGER, INTENT(IN) :: ix
     37  REAL,    INTENT(IN) :: f(:)
     38  REAL,    INTENT(IN) :: x(:)
     39  REAL :: slopes1(SIZE(f,1))
     40!-------------------------------------------------------------------------------
     41! Local:
     42  INTEGER :: n, i, j, sta(2), sto(2)
     43  REAL :: xc(SIZE(f,1))                             ! (n) cell centers
     44  REAL :: h(2:SIZE(f,1)-1), delta_xc(2:SIZE(f,1)-1) ! (2:n-1)
     45  REAL :: fm, ff, fp, dx
     46!-------------------------------------------------------------------------------
     47  n=SIZE(f,ix)
     48  FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2.
     49  FORALL(i=2:n-1)
     50    h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1)
     51  END FORALL
     52  slopes1(:)=0.
     53  DO i=2,n-1
     54    ff=f(i); fm=f(i-1); fp=f(i+1)
     55    IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN
     56      slopes1(i)=0.; CYCLE           !--- Local extremum
     57      !--- 2nd order slope ; (fm, ff, fp) strictly monotonous
     58      slopes1(i)=(fp-fm)/delta_xc(i)
     59      !--- Slope limitation
     60      slopes1(i) = SIGN(MIN(ABS(slopes1(i)), &
     61        ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes1(i) )
     62     END IF
     63  END DO
     64
     65END FUNCTION slopes1
     66!
     67!-------------------------------------------------------------------------------
     68
     69
     70!-------------------------------------------------------------------------------
     71!
     72PURE FUNCTION slopes2(ix, f, x)
     73!
     74!-------------------------------------------------------------------------------
     75! Arguments:
     76  INTEGER, INTENT(IN) :: ix
     77  REAL,    INTENT(IN) :: f(:, :)
     78  REAL,    INTENT(IN) :: x(:)
     79  REAL :: slopes2(SIZE(f,1),SIZE(f,2))
     80!-------------------------------------------------------------------------------
     81! Local:
     82  INTEGER :: n, i, j, sta(2), sto(2)
     83  REAL, ALLOCATABLE :: xc(:)                        ! (n) cell centers
     84  REAL, ALLOCATABLE :: h(:), delta_xc(:)            ! (2:n-1)
     85  REAL :: fm, ff, fp, dx
     86!-------------------------------------------------------------------------------
     87  n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1))
     88  FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2.
     89  FORALL(i=2:n-1)
     90    h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1)
     91  END FORALL
     92  slopes2(:,:)=0.
     93  sta=[1,1]; sta(ix)=2
     94  sto=SHAPE(f); sto(ix)=n-1
     95  DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j)
     96    DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i)
     97      ff=f(i,j)
     98      SELECT CASE(ix)
     99        CASE(1); fm=f(i-1,j); fp=f(i+1,j)
     100        CASE(2); fm=f(i,j-1); fp=f(i,j+1)
     101      END SELECT
     102      IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN
     103        slopes2(i,j)=0.; CYCLE           !--- Local extremum
     104        !--- 2nd order slope ; (fm, ff, fp) strictly monotonous
     105        slopes2(i,j)=(fp-fm)/dx
     106        !--- Slope limitation
     107        slopes2(i,j) = SIGN(MIN(ABS(slopes2(i,j)), &
     108          ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes2(i,j) )
     109       END IF
     110    END DO
     111  END DO
     112  DEALLOCATE(xc,h,delta_xc)
     113
     114END FUNCTION slopes2
     115!
     116!-------------------------------------------------------------------------------
     117
     118
     119!-------------------------------------------------------------------------------
     120!
     121PURE FUNCTION slopes3(ix, f, x)
     122!
     123!-------------------------------------------------------------------------------
     124! Arguments:
     125  INTEGER, INTENT(IN) :: ix
     126  REAL,    INTENT(IN) :: f(:, :, :)
     127  REAL,    INTENT(IN) :: x(:)
     128  REAL :: slopes3(SIZE(f,1),SIZE(f,2),SIZE(f,3))
     129!-------------------------------------------------------------------------------
     130! Local:
     131  INTEGER :: n, i, j, k, sta(3), sto(3)
     132  REAL, ALLOCATABLE :: xc(:)                        ! (n) cell centers
     133  REAL, ALLOCATABLE :: h(:), delta_xc(:)            ! (2:n-1)
     134  REAL :: fm, ff, fp, dx
     135!-------------------------------------------------------------------------------
     136  n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1))
     137  FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2.
     138  FORALL(i=2:n-1)
     139    h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1)
     140  END FORALL
     141  slopes3(:,:,:)=0.
     142  sta=[1,1,1]; sta(ix)=2
     143  sto=SHAPE(f); sto(ix)=n-1
     144  DO k=sta(3),sto(3); IF(ix==3) dx=delta_xc(k)
     145    DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j)
     146      DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i)
     147        ff=f(i,j,k)
     148        SELECT CASE(ix)
     149          CASE(1); fm=f(i-1,j,k); fp=f(i+1,j,k)
     150          CASE(2); fm=f(i,j-1,k); fp=f(i,j+1,k)
     151          CASE(3); fm=f(i,j,k-1); fp=f(i,j,k+1)
     152        END SELECT
     153        IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN
     154          slopes3(i,j,k)=0.; CYCLE           !--- Local extremum
     155          !--- 2nd order slope ; (fm, ff, fp) strictly monotonous
     156          slopes3(i,j,k)=(fp-fm)/dx
     157          !--- Slope limitation
     158          slopes3(i,j,k) = SIGN(MIN(ABS(slopes3(i,j,k)), &
     159            ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes3(i,j,k) )
     160         END IF
     161      END DO
     162    END DO
     163  END DO
     164  DEALLOCATE(xc,h,delta_xc)
     165
     166END FUNCTION slopes3
     167!
     168!-------------------------------------------------------------------------------
     169
     170
     171!-------------------------------------------------------------------------------
     172!
     173PURE FUNCTION slopes4(ix, f, x)
     174!
     175!-------------------------------------------------------------------------------
     176! Arguments:
     177  INTEGER, INTENT(IN) :: ix
     178  REAL,    INTENT(IN) :: f(:, :, :, :)
     179  REAL,    INTENT(IN) :: x(:)
     180  REAL :: slopes4(SIZE(f,1),SIZE(f,2),SIZE(f,3),SIZE(f,4))
     181!-------------------------------------------------------------------------------
     182! Local:
     183  INTEGER :: n, i, j, k, l, m, sta(4), sto(4)
     184  REAL, ALLOCATABLE :: xc(:)                        ! (n) cell centers
     185  REAL, ALLOCATABLE :: h(:), delta_xc(:)            ! (2:n-1)
     186  REAL :: fm, ff, fp, dx
     187!-------------------------------------------------------------------------------
     188  n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1))
     189  FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2.
     190  FORALL(i=2:n-1)
     191    h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1)
     192  END FORALL
     193  slopes4(:,:,:,:)=0.
     194  sta=[1,1,1,1]; sta(ix)=2
     195  sto=SHAPE(f); sto(ix)=n-1
     196  DO l=sta(4),sto(4); IF(ix==4) dx=delta_xc(l)
     197    DO k=sta(3),sto(3); IF(ix==3) dx=delta_xc(k)
     198      DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j)
     199        DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i)
     200          ff=f(i,j,k,l)
     201          SELECT CASE(ix)
     202            CASE(1); fm=f(i-1,j,k,l); fp=f(i+1,j,k,l)
     203            CASE(2); fm=f(i,j-1,k,l); fp=f(i,j+1,k,l)
     204            CASE(3); fm=f(i,j,k-1,l); fp=f(i,j,k+1,l)
     205            CASE(4); fm=f(i,j,k,l-1); fp=f(i,j,k,l+1)
     206          END SELECT
     207          IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN
     208            slopes4(i,j,k,l)=0.; CYCLE           !--- Local extremum
     209            !--- 2nd order slope ; (fm, ff, fp) strictly monotonous
     210            slopes4(i,j,k,l)=(fp-fm)/dx
     211            !--- Slope limitation
     212            slopes4(i,j,k,l) = SIGN(MIN(ABS(slopes4(i,j,k,l)), &
     213              ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes4(i,j,k,l) )
     214           END IF
     215        END DO
     216      END DO
     217    END DO
     218  END DO
     219  DEALLOCATE(xc,h,delta_xc)
     220
     221END FUNCTION slopes4
     222!
     223!-------------------------------------------------------------------------------
     224
     225
     226!-------------------------------------------------------------------------------
     227!
     228PURE FUNCTION slopes5(ix, f, x)
     229!
     230!-------------------------------------------------------------------------------
     231! Arguments:
     232  INTEGER, INTENT(IN) :: ix
     233  REAL,    INTENT(IN) :: f(:, :, :, :, :)
     234  REAL,    INTENT(IN) :: x(:)
     235  REAL :: slopes5(SIZE(f,1),SIZE(f,2),SIZE(f,3),SIZE(f,4),SIZE(f,5))
     236!-------------------------------------------------------------------------------
     237! Local:
     238  INTEGER :: n, i, j, k, l, m, sta(5), sto(5)
     239  REAL, ALLOCATABLE :: xc(:)                        ! (n) cell centers
     240  REAL, ALLOCATABLE :: h(:), delta_xc(:)            ! (2:n-1)
     241  REAL :: fm, ff, fp, dx
     242!-------------------------------------------------------------------------------
     243  n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1))
     244  FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2.
     245  FORALL(i=2:n-1)
     246    h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1)
     247  END FORALL
     248  slopes5(:,:,:,:,:)=0.
     249  sta=[1,1,1,1,1]; sta(ix)=2
     250  sto=SHAPE(f);    sto(ix)=n-1
     251  DO m=sta(5),sto(5); IF(ix==5) dx=delta_xc(m)
     252    DO l=sta(4),sto(4); IF(ix==4) dx=delta_xc(l)
     253      DO k=sta(3),sto(3); IF(ix==3) dx=delta_xc(k)
     254        DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j)
     255          DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i)
     256            ff=f(i,j,k,l,m)
     257            SELECT CASE(ix)
     258              CASE(1); fm=f(i-1,j,k,l,m); fp=f(i+1,j,k,l,m)
     259              CASE(2); fm=f(i,j-1,k,l,m); fp=f(i,j+1,k,l,m)
     260              CASE(3); fm=f(i,j,k-1,l,m); fp=f(i,j,k+1,l,m)
     261              CASE(4); fm=f(i,j,k,l-1,m); fp=f(i,j,k,l+1,m)
     262              CASE(5); fm=f(i,j,k,l,m-1); fp=f(i,j,k,l,m+1)
     263            END SELECT
     264            IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN
     265              slopes5(i,j,k,l,m)=0.; CYCLE           !--- Local extremum
     266              !--- 2nd order slope ; (fm, ff, fp) strictly monotonous
     267              slopes5(i,j,k,l,m)=(fp-fm)/dx
     268              !--- Slope limitation
     269              slopes5(i,j,k,l,m) = SIGN(MIN(ABS(slopes5(i,j,k,l,m)), &
     270                ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes5(i,j,k,l,m) )
     271            END IF
     272          END DO
     273        END DO
     274      END DO
     275    END DO
     276  END DO
     277  DEALLOCATE(xc,h,delta_xc)
     278
     279END FUNCTION slopes5
     280!
     281!-------------------------------------------------------------------------------
     282
     283END MODULE slopes_m

  • LMDZ5/trunk/libf/phylmd/clesphys.h

    r2730 r2788  
    8787       LOGICAL :: ok_qch4
    8888       LOGICAL :: ok_conserv_q
     89       LOGICAL :: adjust_tropopause
     90       LOGICAL :: ok_daily_climoz
    8991
    9092       COMMON/clesphys/                                                 &
    … …  
    130132     &     , iflag_rrtm, ok_strato,ok_hines, ok_qch4                    &
    131133     &     , iflag_ice_thermo, ok_gwd_rando, NSW, iflag_albedo          &
    132      &     , ok_chlorophyll,ok_conserv_q, ok_all_xml
     134     &     , ok_chlorophyll,ok_conserv_q, adjust_tropopause             &
     135     &     , ok_daily_climoz, ok_all_xml
    133136     
    134137       save /clesphys/

  • LMDZ5/trunk/libf/phylmd/conf_phys_m.F90

    r2785 r2788  
    222222    LOGICAL,SAVE  :: carbon_cycle_tr_omp
    223223    LOGICAL,SAVE  :: carbon_cycle_cpl_omp
     224    LOGICAL,SAVE  :: adjust_tropopause_omp
     225    LOGICAL,SAVE  :: ok_daily_climoz_omp
    224226
    225227    INTEGER, INTENT(OUT):: read_climoz ! read ozone climatology, OpenMP shared
    … …  
    20352037    !Config Help = ...
    20362038    !
     2039    !
     2040    adjust_tropopause = .FALSE.
     2041    CALL getin('adjust_tropopause', adjust_tropopause_omp)
     2042    !
     2043    !Config Key  = adjust_tropopause
     2044    !Config Desc = Adjust the ozone field from the climoz file by stretching its
     2045    !              tropopause so that it matches the one of LMDZ.
     2046    !Config Def  = .FALSE.
     2047    !Config Help = Ensure tropospheric ozone column conservation.
     2048    !
     2049    !
     2050    ok_daily_climoz = .TRUE.
     2051    CALL getin('ok_daily_climoz', ok_daily_climoz_omp)
     2052    !
     2053    !Config Key  = ok_daily_climoz
     2054    !Config Desc = Interpolate in time the ozone forcings within ce0l.
     2055    !              .TRUE. if backward compatibility is needed.
     2056    !Config Def  = .TRUE.
     2057    !Config Help = .FALSE. ensure much fewer (no calendar dependency)
     2058    !  and lighter monthly climoz files, inetrpolated in time at gcm run time.
    20372059    !
    20382060    ecrit_LES_omp = 1./8.
    … …  
    22912313    callstats = callstats_omp
    22922314    ecrit_LES = ecrit_LES_omp
     2315    adjust_tropopause = adjust_tropopause_omp
     2316    ok_daily_climoz = ok_daily_climoz_omp
    22932317    carbon_cycle_tr = carbon_cycle_tr_omp
    22942318    carbon_cycle_cpl = carbon_cycle_cpl_omp
    … …  
    24852509    write(lunout,*)' aerosol_couple = ', aerosol_couple
    24862510    write(lunout,*)' flag_aerosol = ', flag_aerosol
    2487     write(lunout,*)' flag_aerosol_strat = ', flag_aerosol_strat
     2511    write(lunout,*)' flag_aerosol_strat= ', flag_aerosol_strat
    24882512    write(lunout,*)' new_aod = ', new_aod
    24892513    write(lunout,*)' aer_type = ',aer_type
    … …  
    25682592    write(lunout,*) 'SSO gkwake=',gkwake
    25692593    write(lunout,*) 'SSO gklift=',gklift
     2594    write(lunout,*) 'adjust_tropopause = ', adjust_tropopause
     2595    write(lunout,*) 'ok_daily_climoz = ',ok_daily_climoz
    25702596    write(lunout,*) 'read_climoz = ', read_climoz
    25712597    write(lunout,*) 'carbon_cycle_tr = ', carbon_cycle_tr

  • LMDZ5/trunk/libf/phylmd/limit_read_mod.F90

    r2770 r2788  
    223223          ierr=NF90_INQUIRE(nid, UnlimitedDimID=ndimid)
    224224          ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
    225           WRITE(abort_message,'(a,2(i0,a))')'limit.nc records number (',nn,') does no'//&
     225          WRITE(abort_message,'(a,2(i3,a))')'limit.nc records number (',nn,') does no'//&
    226226            't match year length (',year_len,')'
    227227          IF(nn/=year_len) CALL abort_physic(modname,abort_message,1)

  • LMDZ5/trunk/libf/phylmd/open_climoz_m.F90

    r1907 r2788  
    66contains
    77
    8   subroutine open_climoz(ncid, press_in_edg)
     8  subroutine open_climoz(ncid, press_in_cen, press_in_edg, time_in)
    99
    1010    ! This procedure should be called once per "gcm" run, by a single
    1111    ! thread of each MPI process.
    1212    ! The root MPI process opens "climoz_LMDZ.nc", reads the pressure
    13     ! levels and broadcasts them to the other processes.
     13    ! levels and the times and broadcasts them to the other processes.
    1414
    1515    ! We assume that, in "climoz_LMDZ.nc", the pressure levels are in hPa
    … …  
    2121    use mod_phys_lmdz_mpi_data, only: is_mpi_root
    2222    use mod_phys_lmdz_mpi_transfert, only: bcast_mpi ! broadcast
     23    USE phys_cal_mod, ONLY: calend, year_len, year_cur
    2324
     25    ! OpenMP shares arguments:
    2426    integer, intent(out):: ncid ! of "climoz_LMDZ.nc", OpenMP shared
    2527
    26     real, pointer:: press_in_edg(:)
    27     ! edges of pressure intervals for ozone climatology, in Pa, in strictly
    28     ! ascending order, OpenMP shared
     28    ! pressure levels for ozone climatology, in Pa, strictly ascending order
     29    real, pointer:: press_in_cen(:) !--- at cells centers
     30    real, pointer:: press_in_edg(:) !--- at the interfaces (pressure intervals)
     31
     32    real, pointer:: time_in(:)
     33    ! times of ozone climatology records, in days since Jan. 1st 0h
    2934
    3035    ! Variables local to the procedure:
    3136
    32     real, pointer:: plev(:)
    33     ! (pressure levels for ozone climatology, converted to Pa, in strictly
    34     ! ascending order)
    35 
    3637    integer varid ! for NetCDF
    3738    integer n_plev ! number of pressure levels in the input data
     39    integer n_time ! number of time records    in the input field
    3840    integer k
     41    CHARACTER(LEN=80)  :: modname
     42    CHARACTER(LEN=320) :: msg
    3943
    4044    !---------------------------------------
    4145
    42     print *, "Call sequence information: open_climoz"
     46    modname="open_climoz"
     47    print *, "Call sequence information: "//TRIM(modname)
    4348
    4449    if (is_mpi_root) then
    … …  
    4651
    4752       call nf95_inq_varid(ncid, "plev", varid)
    48        call nf95_gw_var(ncid, varid, plev)
     53       call nf95_gw_var(ncid, varid, press_in_cen)
    4954       ! Convert from hPa to Pa because "paprs" and "pplay" are in Pa:
    50        plev = plev * 100.
    51        n_plev = size(plev)
     55       press_in_cen = press_in_cen * 100.
     56       n_plev = size(press_in_cen)
     57       CALL NF95_INQ_VARID(ncID, "time", varID)
     58       CALL NF95_GW_VAR(ncid, varid, time_in)
     59       n_time = SIZE(time_in)
     60       IF(ALL([year_len,14]/=n_time)) THEN             !--- Check records number
     61         WRITE(msg,'(a,3(i4,a))')TRIM(modname)//': expected input file could be&
     62         & monthly or daily (14 or ',year_len,' records for year ',year_cur,' a&
     63         &nd calendar "'//TRIM(calend)//'"), but ',n_time,' records were found'
     64         CALL abort_physic(modname, msg, 1)
     65       END IF
     66
    5267    end if
    5368
    5469    call bcast_mpi(n_plev)
    55     if (.not. is_mpi_root) allocate(plev(n_plev))
    56     call bcast_mpi(plev)
    57    
     70    if (.not. is_mpi_root) allocate(press_in_cen(n_plev))
     71    call bcast_mpi(press_in_cen)
     72    call bcast_mpi(n_time)
     73    if (.not. is_mpi_root) allocate(time_in(n_time))
     74    call bcast_mpi(time_in)
     75
    5876    ! Compute edges of pressure intervals:
    5977    allocate(press_in_edg(n_plev + 1))
    … …  
    6179       press_in_edg(1) = 0.
    6280       ! We choose edges halfway in logarithm:
    63        forall (k = 2:n_plev) press_in_edg(k) = sqrt(plev(k - 1) * plev(k))
     81       forall (k = 2:n_plev) press_in_edg(k) = sqrt(press_in_cen(k - 1) * press_in_cen(k))
    6482       press_in_edg(n_plev + 1) = huge(0.)
    6583       ! (infinity, but any value guaranteed to be greater than the
    … …  
    6785    end if
    6886    call bcast_mpi(press_in_edg)
    69     deallocate(plev) ! pointer
     87!    deallocate(press_in_cen) ! pointer
    7088
    7189  end subroutine open_climoz

  • LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90

    r2752 r2788  
    487487!$OMP THREADPRIVATE(budg_sed_part)
    488488#endif
     489      REAL, ALLOCATABLE, SAVE, DIMENSION(:) :: pr_tropopause
     490!$OMP THREADPRIVATE(pr_tropopause)
    489491
    490492CONTAINS
    … …  
    764766      ALLOCATE (vsed_aer(klon,klev))
    765767#endif
     768      ALLOCATE (pr_tropopause(klon))
    766769
    767770END SUBROUTINE phys_local_var_init
    … …  
    10181021      DEALLOCATE (budg_sed_part)
    10191022#endif
     1023      DEALLOCATE (pr_tropopause)
    10201024
    10211025END SUBROUTINE phys_local_var_end

  • LMDZ5/trunk/libf/phylmd/physiq_mod.F90

    r2784 r2788  
    3535    USE change_srf_frac_mod
    3636    USE surface_data,     ONLY : type_ocean, ok_veget, ok_snow
     37    USE tropopause_m,     ONLY: dyn_tropopause
    3738#ifdef CPP_Dust
    3839    USE phytracr_spl_mod, ONLY: phytracr_spl
    … …  
    191192       beta_prec,  &
    192193       rneb,  &
    193        zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic
     194       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic, &
     195       pr_tropopause
    194196       !
    195197    USE phys_state_var_mod ! Variables sauvegardees de la physique
    … …  
    205207    USE phys_output_ctrlout_mod
    206208    use open_climoz_m, only: open_climoz ! ozone climatology from a file
    207     use regr_pr_av_m, only: regr_pr_av
     209    use regr_pr_time_av_m, only: regr_pr_time_av
    208210    use netcdf95, only: nf95_close
    209211    !IM for NMC files
    … …  
    10061008    !$OMP THREADPRIVATE(first)
    10071009
    1008     integer, save::  read_climoz ! read ozone climatology
     1010    ! VARIABLES RELATED TO OZONE CLIMATOLOGIES ; all are OpenMP shared
     1011    ! Note that pressure vectors are in Pa and in stricly ascending order
     1012    INTEGER,SAVE :: read_climoz                ! Read ozone climatology
    10091013    !     (let it keep the default OpenMP shared attribute)
    10101014    !     Allowed values are 0, 1 and 2
    … …  
    10131017    !     2: read two ozone climatologies, the average day and night
    10141018    !     climatology and the daylight climatology
    1015 
    1016     integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
    1017     !     (let it keep the default OpenMP shared attribute)
    1018 
    1019     real, pointer, save:: press_climoz(:)
    1020     ! (let it keep the default OpenMP shared attribute)
    1021     ! edges of pressure intervals for ozone climatologies, in Pa, in strictly
    1022     ! ascending order
    1023 
    1024     integer ro3i
    1025     !     required time index in NetCDF file for the ozone fields, between 1
    1026     !     and 360
     1019    INTEGER,SAVE :: ncid_climoz                ! NetCDF file identifier
     1020    REAL, POINTER, SAVE :: press_cen_climoz(:) ! Pressure levels
     1021    REAL, POINTER, SAVE :: press_edg_climoz(:) ! Edges of pressure intervals
     1022    REAL, POINTER, SAVE :: time_climoz(:)      ! Time vector
     1023    CHARACTER(LEN=13), PARAMETER :: vars_climoz(2) &
     1024                                  = ["tro3         ","tro3_daylight"]
     1025    ! vars_climoz(1:read_climoz): variables names in climoz file.
    10271026
    10281027    INTEGER ierr
    … …  
    16711670
    16721671       !$omp single
    1673        IF (read_climoz >= 1) THEN
    1674           CALL open_climoz(ncid_climoz, press_climoz)
    1675        ENDIF
     1672       IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz,   &
     1673                                         press_edg_climoz, time_climoz)
    16761674       !$omp end single
    16771675       !
    … …  
    19781976          wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
    19791977       ELSE
    1980           ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1
    1981 
    1982           IF (ro3i == 361) ro3i = 360
    1983           ! (This should never occur, except perhaps because of roundup
    1984           ! error. See documentation.)
    1985 
    1986           IF (read_climoz == 1) THEN
    1987              CALL regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
    1988                   press_in_edg=press_climoz, paprs=paprs, v3=wo)
     1978          IF(adjust_tropopause) THEN
     1979            WRITE(*,*)' >> Adjusting O3 field to match gcm tropopause.'
     1980            CALL dyn_tropopause(t_seri, ztsol, paprs, pplay, presnivs, rot,  &
     1981                 pr_tropopause)
     1982            CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),    &
     1983                 days_elapsed+jh_cur-jh_1jan, press_edg_climoz, paprs, wo,   &
     1984                 time_climoz,  latitude_deg,  press_cen_climoz, pr_tropopause)
    19891985          ELSE
    1990              ! read_climoz == 2
    1991              CALL regr_pr_av(ncid_climoz, (/"tro3         ", &
    1992                   "tro3_daylight"/), julien=ro3i, press_in_edg=press_climoz, &
    1993                   paprs=paprs, v3=wo)
    1994           ENDIF
     1986            WRITE(*,*)' >> Interpolating O3 field directly on gcm levels.'
     1987            CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),    &
     1988                 days_elapsed+jh_cur-jh_1jan, press_edg_climoz, paprs, wo,  &
     1989                 time_climoz)
     1990          END IF
    19951991          ! Convert from mole fraction of ozone to column density of ozone in a
    19961992          ! cell, in kDU:
    19971993          FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
    19981994               * zmasse / dobson_u / 1e3
    1999           ! (By regridding ozone values for LMDZ only once per day, we
     1995          ! (By regridding ozone values for LMDZ only once a day, we
    20001996          ! have already neglected the variation of pressure in one
    20011997          ! day. So do not recompute "wo" at each time step even if
    … …  
    44654461             CALL nf95_close(ncid_climoz)
    44664462          ENDIF
    4467           DEALLOCATE(press_climoz) ! pointer
     4463          DEALLOCATE(press_edg_climoz) ! pointer
     4464          DEALLOCATE(press_cen_climoz) ! pointer
    44684465       ENDIF
    44694466       !$OMP END MASTER

  • LMDZ5/trunk/libf/phylmd/regr_lat_time_coefoz_m.F90

    r2440 r2788  
    4141
    4242    use mod_grid_phy_lmdz, ONLY : nbp_lat
    43     use regr1_conserv_m, only: regr1_conserv
    44     use regr3_lint_m, only: regr3_lint
     43    use regr_conserv_m, only: regr_conserv
     44    use regr_lint_m, only: regr_lint
    4545    use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, handle_err, &
    4646         nf95_put_var, nf95_gw_var
    … …  
    8484    ! Last dimension is month number.)
    8585
    86     real, allocatable:: o3_par_out(:, :, :) ! (jjm + 1, n_plev, 360)
     86    real, allocatable:: o3_par_out(:, :, :) ! (jjm + 1, n_plev, ndays)
    8787    ! (regridded ozone parameter)
    8888    ! ("o3_par_out(j, l, day)" is at latitude "rlatu(j)", pressure
    … …  
    162162    latitude = latitude / 180. * pi
    163163    n_lat = size(latitude)
    164     ! We need to supply the latitudes to "regr1_conserv" in
     164    ! We need to supply the latitudes to "regr_conserv" in
    165165    ! ascending order, so invert order if necessary:
    166166    desc_lat = latitude(1) > latitude(n_lat)
    … …  
    209209       ! We average with respect to sine of latitude, which is
    210210       ! equivalent to weighting by cosine of latitude:
    211        call regr1_conserv(o3_par_in, xs = sin(lat_in_edg), &
     211       call regr_conserv(1, o3_par_in, xs = sin(lat_in_edg), &
    212212            xt = (/-1., sin((/boundslat_reg(nbp_lat-1:1:-1,south)/)), 1./), &
    213213            vt = v_regr_lat(nbp_lat:1:-1, :, 1:12))
    … …  
    221221
    222222       ! Regrid in time by linear interpolation:
    223        o3_par_out = regr3_lint(v_regr_lat, tmidmonth, tmidday)
     223       call regr_lint(3, v_regr_lat, tmidmonth, tmidday, o3_par_out)
    224224
    225225       ! Write to file:

  • LMDZ5/trunk/libf/phylmd/regr_pr_comb_coefoz_m.F90

    r2346 r2788  
    7676    use dimphy, only: klon
    7777    use mod_phys_lmdz_mpi_data, only: is_mpi_root
    78     use regr_pr_av_m, only: regr_pr_av
     78    use regr_pr_time_av_m, only: regr_pr_time_av
    7979    use regr_pr_int_m, only: regr_pr_int
    8080    use press_coefoz_m, only: press_in_edg, plev
    … …  
    128128    !$omp end master
    129129
    130     call regr_pr_av(ncid, (/"a2       ", "a4       ", "a6       ", &
    131          "P_net_Mob", "r_Mob    ", "temp_Mob ", "R_Het    "/), julien, &
     130    call regr_pr_time_av(ncid, (/"a2       ", "a4       ", "a6       ", &
     131         "P_net_Mob", "r_Mob    ", "temp_Mob ", "R_Het    "/), REAL(julien-1), &
    132132         press_in_edg, paprs, coefoz)
    133133    a2 = coefoz(:, :, 1)

  • LMDZ5/trunk/libf/phylmd/regr_pr_int_m.F90

    r2346 r2788  
    2828    use netcdf, only: nf90_get_var
    2929    use assert_m, only: assert
    30     use regr1_lint_m, only: regr1_lint
     30    use regr_lint_m, only: regr_lint
    3131    use mod_phys_lmdz_mpi_data, only: is_mpi_root
    3232    use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, nbp_lev
    … …  
    9797    ! Regrid in pressure at each horizontal position:
    9898    do i = 1, klon
    99        v3(i, nbp_lev:1:-1) = regr1_lint(v2(i, :), (/0., plev/), pplay(i, nbp_lev:1:-1))
     99       call regr_lint(1, v2(i, :), (/0., plev/), pplay(i, nbp_lev:1:-1), &
     100                         v3(i, nbp_lev:1:-1))
    100101       ! (invert order of indices because "pplay" is in descending order)
    101102    end do

  • LMDZ5/trunk/libf/phylmd/regr_pr_o3_m.F90

    r2440 r2788  
    2828    use netcdf, only:  nf90_nowrite, nf90_get_var
    2929    use assert_m, only: assert
    30     use regr1_conserv_m, only: regr1_conserv
     30    use regr_conserv_m, only: regr_conserv
    3131    use press_coefoz_m, only: press_in_edg
    3232    use time_phylmdz_mod, only: day_ref
    … …  
    7575    ! Poles:
    7676    do j = 1, nbp_lat, nbp_lat-1
    77        call regr1_conserv(r_mob(j, :), press_in_edg, &
     77       call regr_conserv(1, r_mob(j, :), press_in_edg, &
    7878            p3d(1, j, nbp_lev + 1:1:-1), o3_mob_regr(1, j, nbp_lev:1:-1))
    7979       ! (invert order of indices because "p3d" is in descending order)
    … …  
    8383    do j = 2, nbp_lat-1
    8484       do i = 1, nbp_lon
    85           call regr1_conserv(r_mob(j, :), press_in_edg, &
     85          call regr_conserv(1, r_mob(j, :), press_in_edg, &
    8686               p3d(i, j, nbp_lev + 1:1:-1), o3_mob_regr(i, j, nbp_lev:1:-1))
    8787          ! (invert order of indices because "p3d" is in descending order)
Note: See TracChangeset for help on using the changeset viewer.