Ignore:
Timestamp:
Mar 30, 2017, 4:16:38 PM (8 years ago)
Author:
Laurent Fairhead
Message:

Merged trunk changes r2785:2838 into testing branch

Location:
LMDZ5/branches/testing
Files:
9 deleted
52 edited
12 copied

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/libf/dyn3d_common/disvert.F90

    r2641 r2839  
    1010  use new_unit_m, only: new_unit
    1111  use assert_m, only: assert
    12   USE comvert_mod, ONLY: ap, bp, nivsigs, nivsig, dpres, presnivs, &
    13                          pa, preff, scaleheight
     12  USE comvert_mod, ONLY: ap, bp, aps, bps, nivsigs, nivsig, dpres, presnivs, &
     13                         pseudoalt, pa, preff, scaleheight
    1414  USE logic_mod, ONLY: ok_strato
    1515
    … …  
    346346  DO l = 1, llm
    347347     dpres(l) = bp(l) - bp(l+1)
     348     aps(l) =  0.5 *( ap(l) +ap(l+1))
     349     bps(l) =  0.5 *( bp(l) +bp(l+1))
    348350     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
     351     pseudoalt(l) = log(preff/presnivs(l))*scaleheight
    349352     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
    350           log(preff/presnivs(l))*scaleheight &
     353          pseudoalt(l) &
    351354          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
    352355          max(ap(l+1)+bp(l+1)*preff, 1.e-10))

  • LMDZ5/branches/testing/libf/dyn3dmem/dynetat0_loc.f90

    r2641 r2839  
    153153  DO iq=1,nqtot
    154154    var=tname(iq)
     155#ifdef INCA
     156    IF (var .eq. "O3" ) THEN
     157       IF(NF90_INQ_VARID(fID,var,vID) == NF90_NoErr) THEN
     158          CALL get_var2(var,q_glo); q(ijb_u:ije_u,:,iq)=q_glo(ijb_u:ije_u,:); CYCLE
     159       ELSE
     160          WRITE(lunout,*) 'Tracer O3 is missing - it is initialized to OX'
     161          IF(NF90_INQ_VARID(fID,"OX",vID) == NF90_NoErr) THEN
     162             CALL get_var2("OX",q_glo); q(ijb_u:ije_u,:,iq)=q_glo(ijb_u:ije_u,:); CYCLE
     163          ENDIF
     164       ENDIF
     165    ENDIF
     166#endif
    155167    IF(NF90_INQ_VARID(fID,var,vID)==NF90_NoErr) THEN
    156168      CALL get_var2(var,q_glo); q(ijb_u:ije_u,:,iq)=q_glo(ijb_u:ije_u,:); CYCLE

  • LMDZ5/branches/testing/libf/dynphy_lonlat/inigeomphy_mod.F90

    r2594 r2839  
    2424  USE mod_interface_dyn_phys, ONLY :  init_interface_dyn_phys
    2525  USE nrtype, ONLY: pi
     26  USE comvert_mod, ONLY: preff, ap, bp, aps, bps, presnivs, &
     27                         scaleheight, pseudoalt
     28  USE vertical_layers_mod, ONLY: init_vertical_layers
    2629  IMPLICIT NONE
    2730
    … …  
    216219                     airefi,cufi,cvfi)
    217220
     221  ! copy over preff , ap(), bp(), etc
     222  CALL init_vertical_layers(nlayer,preff,scaleheight, &
     223                            ap,bp,aps,bps,presnivs,pseudoalt)
     224
    218225!$OMP END PARALLEL
    219226

  • LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90

    r2787 r2839  
    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/branches/testing/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90

    r2669 r2839  
    1212                     prad,pg,pr,pcpp,iflag_phys)
    1313  USE dimphy, ONLY: init_dimphy
    14   USE comvert_mod, ONLY: preff, ap, bp, presnivs, scaleheight, pseudoalt
    1514  USE inigeomphy_mod, ONLY: inigeomphy
    1615  USE mod_grid_phy_lmdz, ONLY: nbp_lon,nbp_lat,nbp_lev,klon_glo ! number of atmospheric columns (on full grid)
    … …  
    110109!$OMP COPYIN(annee_ref, day_ini, day_ref, start_time)
    111110
    112   ! copy over preff , ap(), bp(), etc
    113   CALL init_vertical_layers(nlayer,preff,scaleheight, &
    114                             ap,bp,presnivs,pseudoalt)
    115 
    116111  ! Initialize physical constants in physics:
    117112  CALL inifis(punjours,prad,pg,pr,pcpp)

  • LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/init_ssrf_m.F90

    r2641 r2839  
    99  USE grid_atob_m,        ONLY: grille_m
    1010  USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo
     11  USE ioipsl_getin_p_mod, ONLY: getin_p
    1112  USE comconst_mod, ONLY: im, pi
    1213
    … …  
    4243  REAL, ALLOCATABLE :: dlat_lic(:), lat_lic(:,:), flic_tmp(:,:), vtmp(:,:)
    4344  REAL              :: date, lev(1), dt, deg2rad
     45  LOGICAL           :: no_ter_antartique   ! If true, no land points are allowed at Antartic
    4446!-------------------------------------------------------------------------------
    4547  deg2rad= pi/180.0
    … …  
    9799  END DO
    98100
     101
     102  !--- Option no_ter_antartique removes all land fractions souther than 60S.
     103  !--- Land ice is set instead of the land fractions on these latitudes.
     104  !--- The ocean and sea-ice fractions are not changed.
     105  no_ter_antartique=.FALSE.
     106  CALL getin_p('no_ter_antartique',no_ter_antartique)
     107  WRITE(lunout,*)"no_ter_antartique=",no_ter_antartique
     108  IF (no_ter_antartique) THEN
     109     ! Remove all land fractions souther than 60S and set land-ice instead
     110     WRITE(lunout,*) "Remove land fractions souther than 60deg south by increasing"
     111     WRITE(lunout,*) "the continental ice fractions. No land can now be found at Antartic."
     112     DO ji=1, klon
     113        IF (latitude_deg(ji)<-60.0) THEN
     114           pctsrf(ji,is_lic) = pctsrf(ji,is_lic) + pctsrf(ji,is_ter)
     115           pctsrf(ji,is_ter) = 0
     116        END IF
     117     END DO
     118  END IF
     119
     120
    99121!--- Sub-surface ocean and sea ice (sea ice set to zero for start).
    100122  pctsrf(:,is_oce)=(1.-zmasq(:))

  • LMDZ5/branches/testing/libf/dynphy_lonlat/phymar/iniphysiq_mod.F90

    r2641 r2839  
    1212                     prad,pg,pr,pcpp,iflag_phys)
    1313  USE dimphy, ONLY: init_dimphy
    14   USE comvert_mod, ONLY: preff, ap, bp, presnivs, scaleheight, pseudoalt
    1514  USE inigeomphy_mod, ONLY: inigeomphy
    1615  USE vertical_layers_mod, ONLY : init_vertical_layers
    … …  
    7170!$OMP PARALLEL
    7271
    73   ! copy over preff , ap(), bp(), etc
    74   CALL init_vertical_layers(nlayer,preff,scaleheight, &
    75                             ap,bp,presnivs,pseudoalt)
    76 
    7772  ! Initialize tracer names, numbers, etc. for physics
    7873  CALL init_infotrac_phy(nqtot)

  • LMDZ5/branches/testing/libf/misc/slopes_m.F90

    r2471 r2839  
    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/branches/testing/libf/phylmd/YOETHF.h

    r2641 r2839  
    1919      REAL R5ALVCP,R5ALSCP,RALVDCP,RALSDCP,RALFDCP,RTWAT,RTBER,RTBERCU
    2020      REAL RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,RKOOP2
     21      LOGICAL OK_BAD_ECMWF_THERMO ! If TRUE, then variables set by rrtm/suphec.F90
     22                                  ! If FALSE, then variables set by suphel.F90
    2123      COMMON /YOETHF/R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES,    &
    2224     &               RVTMP2, RHOH2O,                                    &
    … …  
    2426     &               RALFDCP,RTWAT,RTBER,RTBERCU,                       &
    2527     &               RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,&
    26      &               RKOOP2
     28     &               RKOOP2,                                            &
     29     &               OK_BAD_ECMWF_THERMO
    2730
    2831!$OMP THREADPRIVATE(/YOETHF/)

  • LMDZ5/branches/testing/libf/phylmd/aero_mod.F90

    r2594 r2839  
    4646! 2/ Total number of aerosols for which an aerosol mass is provided
    4747
    48   INTEGER, PARAMETER :: naero_spc = 10
     48  INTEGER, PARAMETER :: naero_spc = 13
    4949
    5050! Corresponding names for the aerosols
    … …  
    5959       "CIDUSTM", &
    6060       "AIBCM  ", &
    61        "AIPOMM " /)
     61       "AIPOMM " ,&
     62       "ASNO3M ", &
     63       "CSNO3M ", &
     64       "CINO3M " /)
    6265
    6366! 3/ Number of aerosol groups

  • LMDZ5/branches/testing/libf/phylmd/calcul_divers.h

    r2408 r2839  
    55      IF(itap.EQ.1) THEN
    66         itapm1=0
    7 !        surface terre
    8          DO i=1, klon
    9             IF(pctsrf(i,is_ter).GT.0.) THEN
    10                paire_ter(i)=cell_area(i)*pctsrf(i,is_ter)
    11             ENDIF
    12          ENDDO
    137      ENDIF
    148

  • LMDZ5/branches/testing/libf/phylmd/clesphys.h

    r2787 r2839  
    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/branches/testing/libf/phylmd/climb_hq_mod.F90

    r2187 r2839  
    99  SAVE
    1010  PRIVATE
    11   PUBLIC :: climb_hq_down, climb_hq_up
     11  PUBLIC :: climb_hq_down, climb_hq_up, d_h_col_vdf, f_h_bnd
    1212
    1313  REAL, DIMENSION(:,:), ALLOCATABLE :: gamaq, gamah
    … …  
    2323  REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefhq
    2424  !$OMP THREADPRIVATE(Kcoefhq)
     25  REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: h_old ! for diagnostics, h before solving diffusion
     26  !$OMP THREADPRIVATE(h_old)
     27  REAL, SAVE, DIMENSION(:), ALLOCATABLE :: d_h_col_vdf ! for diagnostics, vertical integral of enthalpy change
     28  !$OMP THREADPRIVATE(d_h_col_vdf)
     29  REAL, SAVE, DIMENSION(:), ALLOCATABLE :: f_h_bnd ! for diagnostics, enthalpy flux at surface
     30  !$OMP THREADPRIVATE(f_h_bnd)
    2531
    2632CONTAINS
    … …  
    7177    LOGICAL, SAVE                            :: first=.TRUE.
    7278    !$OMP THREADPRIVATE(first)
    73     REAL, DIMENSION(klon,klev)               :: local_H
     79! JLD now renamed h_old and declared in module
     80!    REAL, DIMENSION(klon,klev)               :: local_H
    7481    REAL, DIMENSION(klon)                    :: psref
    7582    REAL                                     :: delz, pkh
    7683    INTEGER                                  :: k, i, ierr
    77 
    7884! Include
    7985!****************************************************************************************
    … …  
    113119       ALLOCATE(gamah(1:klon,2:klev), STAT=ierr)
    114120       IF ( ierr /= 0 ) PRINT*,' pb in allloc gamah, ierr=', ierr
     121       
     122       ALLOCATE(h_old(klon,klev), STAT=ierr)
     123       IF ( ierr /= 0 )  PRINT*,' pb in allloc h_old, ierr=', ierr
     124       
     125       ALLOCATE(d_h_col_vdf(klon), STAT=ierr)
     126       IF ( ierr /= 0 )  PRINT*,' pb in allloc d_h_col_vdf, ierr=', ierr
     127       
     128       ALLOCATE(f_h_bnd(klon), STAT=ierr)
     129       IF ( ierr /= 0 )  PRINT*,' pb in allloc f_h_bnd, ierr=', ierr
    115130    END IF
    116131
    … …  
    177192!
    178193!****************************************************************************************
    179     local_H(:,:) = 0.0
     194    h_old(:,:) = 0.0
    180195
    181196    DO k=1,klev
    182197       DO i = 1, knon
    183198          ! convertie la temperature en entalpie potentielle
    184           local_H(i,k) = RCPD * temp(i,k) * &
     199          h_old(i,k) = RCPD * temp(i,k) * &
    185200               (psref(i)/pplay(i,k))**RKAPPA
    186201       ENDDO
    187202    ENDDO
    188203
    189     CALL calc_coef(knon, Kcoefhq(:,:), gamah(:,:), delp(:,:), local_H(:,:), &
     204    CALL calc_coef(knon, Kcoefhq(:,:), gamah(:,:), delp(:,:), h_old(:,:), &
    190205         Ccoef_H(:,:), Dcoef_H(:,:), Acoef_H, Bcoef_H)
    191206 
    … …  
    349364! 1)
    350365! Definition of some variables
     366    REAL, DIMENSION(klon,klev)               :: d_h, zairm
    351367!
    352368!****************************************************************************************
    … …  
    355371    d_q(:,:)    = 0.0
    356372    d_t(:,:)    = 0.0
     373    d_h(:,:)    = 0.0
     374    f_h_bnd(:)= 0.0
    357375
    358376    psref(1:knon) = paprs(1:knon,1) 
    … …  
    393411    q_new(1:knon,1) = Acoef_Q(1:knon) + Bcoef_Q(1:knon)*flx_q1(1:knon)*dtime
    394412    h_new(1:knon,1) = Acoef_H(1:knon) + Bcoef_H(1:knon)*flx_h1(1:knon)*dtime
    395    
     413    f_h_bnd(1:knon) = flx_h1(1:knon)
    396414!- All the other layers
    397415    DO k = 2, klev
    … …  
    427445!
    428446!****************************************************************************************
    429 
     447    d_h_col_vdf(:) = 0.0
    430448    DO k = 1, klev
    431449       DO i = 1, knon
    432450          d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k)
    433451          d_q(i,k) = q_new(i,k) - q_old(i,k)
    434        END DO
     452          d_h(i,k) = h_new(i,k) - h_old(i,k)
     453!JLD          d_t(i,k) = d_h(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD    !correction a venir
     454    ! layer air mass
     455          zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
     456          d_h_col_vdf(i) = d_h_col_vdf(i) + d_h(i,k)*zairm(i,k)
     457        END DO
    435458    END DO
    436459
    … …  
    448471       DEALLOCATE(Kcoefhq,stat=ierr)
    449472       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
     473       DEALLOCATE(h_old, d_h_col_vdf, f_h_bnd, stat=ierr)
     474       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate h_old, d_h_col_vdf, f_h_bnd, ierr=', ierr
    450475    END IF
    451476  END SUBROUTINE climb_hq_up

  • LMDZ5/branches/testing/libf/phylmd/concvl.F90

    r2488 r2839  
    55                  d_t, d_q, d_u, d_v, d_tra, &
    66                  rain, snow, kbas, ktop, sigd, &
    7                   cbmf, plcl, plfc, wbeff, upwd, dnwd, dnwdbis, &
     7                  cbmf, plcl, plfc, wbeff, convoccur, &
     8                  upwd, dnwd, dnwdbis, &
    89                  Ma, mip, Vprecip, &
    910                  cape, cin, tvp, Tconv, iflag, &
    … …  
    157158  REAL qs(klon, klev), qs_wake(klon, klev)
    158159  REAL cbmf(klon), plcl(klon), plfc(klon), wbeff(klon)
     160  REAL convoccur(klon)
    159161!LF          SAVE cbmf
    160162!IM/JYG      REAL, SAVE, ALLOCATABLE :: cbmf(:)
    … …  
    228230    ALLOCATE (t1(klon,klev))
    229231    ALLOCATE (q1(klon,klev))
     232!
     233    convoccur(:) = 0.
     234!
    230235    itap = 0
    231236    igout = klon/2 + 1/klon
    … …  
    450455  END DO
    451456
     457  IF (iflag_con==3) THEN
     458    DO i = 1,klon
     459      IF (wbeff(i) > 100. .OR. wbeff(i) == 0 .OR. iflag(i) > 3) THEN
     460        wbeff(i) = 0.
     461        convoccur(i) = 0. 
     462      ELSE
     463        convoccur(i) = 1.
     464      ENDIF
     465    ENDDO
     466  ENDIF
     467
    452468  IF (iflag_con==30) THEN
    453469    DO itra = 1, ntra

  • LMDZ5/branches/testing/libf/phylmd/conf_phys_m.F90

    r2787 r2839  
    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 = .FALSE.
     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/branches/testing/libf/phylmd/cosp/cosp_output_mod.F90

    r2720 r2839  
    2020      INTEGER, DIMENSION(3), SAVE  :: nhoricosp,nvert,nvertmcosp,nvertcol,nvertbze, &
    2121                                      nvertsratio,nvertisccp,nvertp,nverttemp,nvertmisr, &
    22                                       nvertReffIce,nvertReffLiq
     22                                      nvertReffIce,nvertReffLiq,nverttau
    2323      REAL, DIMENSION(3), SAVE                :: zoutm_cosp
    2424!$OMP THREADPRIVATE(nhoricosp, nvert,nvertmcosp,nvertcol,nvertsratio,nvertbze,nvertisccp,nvertp,zoutm_cosp,nverttemp,nvertmisr)
    25 !$OMP THREADPRIVATE(nvertReffIce,nvertReffLiq)
     25!$OMP THREADPRIVATE(nvertReffIce,nvertReffLiq,nverttau)
    2626      REAL, SAVE                   :: zdtimemoy_cosp
    2727!$OMP THREADPRIVATE(zdtimemoy_cosp)
    … …  
    207207  SUBROUTINE cosp_output_open(Nlevlmdz, Ncolumns, presnivs, dtime, freq_cosp, &
    208208                              ok_mensuelCOSP, ok_journeCOSP, ok_hfCOSP, ok_all_xml,  &
    209                               ecrit_mth, ecrit_day, ecrit_hf, use_vgrid, vgrid)
    210 
     209                              ecrit_mth, ecrit_day, ecrit_hf, use_vgrid, vgrid, stlidar)
    211210
    212211  USE iophy
    … …  
    229228  logical                  :: ok_mensuelCOSP, ok_journeCOSP, ok_hfCOSP, use_vgrid, ok_all_xml                   
    230229  type(cosp_vgrid)   :: vgrid   ! Information on vertical grid of stats
     230  type(cosp_lidarstats) :: stlidar ! Summary statistics from lidar simulator
    231231
    232232!!! Variables locales
    … …  
    236236  real,dimension(2,SR_BINS) :: sratio_bounds
    237237  real,dimension(SR_BINS)   ::  sratio_ax
     238  real,dimension(DBZE_BINS) ::  dbze_ax
    238239  CHARACTER(LEN=20), DIMENSION(3)  :: chfreq = (/ '1day', '1d  ', '3h  ' /)           
    239240
    … …  
    257258    enddo
    258259
    259 !    do ii=1,DBZE_BINS
    260 !     dbze_ax(i) = CFAD_ZE_MIN + CFAD_ZE_WIDTH*(ii - 0.5)
    261 !    enddo
    262 
    263 !   sratio_bounds(2,:)=stlidar%srbval(:) ! srbval contains the upper
     260    do i=1,DBZE_BINS
     261     dbze_ax(i) = CFAD_ZE_MIN + CFAD_ZE_WIDTH*(i - 0.5)
     262    enddo
     263 
     264
     265    sratio_bounds(2,:)=stlidar%srbval(:) ! srbval contains the upper
    264266!                                         limits from lmd_ipsl_stats.f90
    265 !   sratio_bounds(1,2:SR_BINS) = stlidar%srbval(1:SR_BINS-1)
    266 !   sratio_bounds(1,1)         = 0.0
    267 !   sratio_bounds(2,SR_BINS)   = 1.e5 ! This matches with Chepfer et al., JGR,
     267    sratio_bounds(1,2:SR_BINS) = stlidar%srbval(1:SR_BINS-1)
     268    sratio_bounds(1,1)         = 0.0
     269    sratio_bounds(2,SR_BINS)   = 1.e5 ! This matches with Chepfer et al., JGR,
    268270!                                    ! 2009. However, it is not consistent
    269271                                     ! with the upper limit in
    270272                                     ! lmd_ipsl_stats.f90, which is
    271273                                     ! LIDAR_UNDEF-1=998.999
    272 !    sratio_ax(:) = (sratio_bounds(1,:)+sratio_bounds(2,:))/2.0
     274     sratio_ax(:) = (sratio_bounds(1,:)+sratio_bounds(2,:))/2.0
    273275
    274276    cosp_outfilenames(1) = 'histmthCOSP'
    … …  
    331333   CALL wxios_add_vaxis("temp", LIDAR_NTEMP, LIDAR_PHASE_TEMP)
    332334   CALL wxios_add_vaxis("cth16", MISR_N_CTH, MISR_CTH)
    333 !   CALL wxios_add_vaxis("dbze", DBZE_BINS, dbze_ax)
    334 !   CALL wxios_add_vaxis("scatratio", SR_BINS, sratio_ax)
     335   CALL wxios_add_vaxis("dbze", DBZE_BINS, dbze_ax)
     336   CALL wxios_add_vaxis("scatratio", SR_BINS, sratio_ax)
    335337   CALL wxios_add_vaxis("ReffIce", numMODISReffIceBins, reffICE_binCenters)
    336338   CALL wxios_add_vaxis("ReffLiq", numMODISReffLiqBins, reffLIQ_binCenters)
     339   print*,'reffICE_binCenters=',reffICE_binCenters
     340   CALL wxios_add_vaxis("tau", 7, ISCCP_TAU)
    337341
    338342#endif
    … …  
    384388                    nvertReffLiq(iff))
    385389
    386 !      CALL histvert(cosp_nidfiles(iff),"dbze","equivalent_reflectivity_factor","dBZ",DBZE_BINS,dbze_ax,nvertbze(iff))
     390      CALL histvert(cosp_nidfiles(iff),"dbze","equivalent_reflectivity_factor","dBZ",DBZE_BINS,dbze_ax,nvertbze(iff))
    387391     
    388 !      CALL histvert(cosp_nidfiles(iff),"scatratio","backscattering_ratio","1",SR_BINS,sratio_ax,nvertsratio(iff))
     392      CALL histvert(cosp_nidfiles(iff),"scatratio","backscattering_ratio","1",SR_BINS,sratio_ax,nvertsratio(iff))
     393
     394      CALL histvert(cosp_nidfiles(iff),"tau","cloud optical depth","1",7,ISCCP_TAU,nverttau(iff))
    389395     
    390396!!! Valeur indefinie en cas IOIPSL

  • LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_write_mod.F90

    r2787 r2839  
    1717   CONTAINS
    1818
    19   SUBROUTINE cosp_output_write(Nlevlmdz, Npoints, Ncolumns, itap, dtime, freq_COSP, &
     19  SUBROUTINE cosp_output_write(Nlevlmdz, Npoints, Ncolumns, itap, dtime, freq_COSP, missing_cosp, &
    2020                               cfg, gbx, vgrid, sglidar, sgradar, stlidar, stradar, &
    2121                               isccp, misr, modis)
    … …  
    3232!!! Variables d'entree
    3333  integer               :: itap, Nlevlmdz, Ncolumns, Npoints
    34   real                  :: freq_COSP, dtime
     34  real                  :: freq_COSP, dtime, missing_val, missing_cosp
    3535  type(cosp_config)     :: cfg     ! Control outputs
    3636  type(cosp_gridbox)    :: gbx     ! Gridbox information. Input for COSP
    … …  
    4949  integer               :: itau_wcosp
    5050  real, dimension(Npoints,PARASOL_NREFL) :: parasolcrefl, Ncref
     51
     52
     53#ifdef CPP_XIOS
     54  missing_val=missing_cosp
     55#else
     56  missing_val=0.
     57#endif
    5158
    5259  Nlevout = vgrid%Nlvgrid
    … …  
    95102     do ip = 1,Npoints
    96103     if(stlidar%lidarcld(ip,k).eq.R_UNDEF)then
    97       stlidar%lidarcld(ip,k)=Cosp_fill_value
     104      stlidar%lidarcld(ip,k)=missing_val
    98105     endif
    99106     enddo
    … …  
    102109      do ip = 1,Npoints
    103110       if(stlidar%cfad_sr(ip,ii,k).eq.R_UNDEF)then
    104         stlidar%cfad_sr(ip,ii,k)=Cosp_fill_value
     111        stlidar%cfad_sr(ip,ii,k)=missing_val
    105112       endif
    106113      enddo
    … …  
    111118   do k = 1,Nlevlmdz
    112119     if(sglidar%beta_mol(ip,k).eq.R_UNDEF)then
    113       sglidar%beta_mol(ip,k)=Cosp_fill_value
     120      sglidar%beta_mol(ip,k)=missing_val
    114121     endif
    115122
    116123     do ii= 1,Ncolumns
    117124       if(sglidar%beta_tot(ip,ii,k).eq.R_UNDEF)then
    118         sglidar%beta_tot(ip,ii,k)=Cosp_fill_value
     125        sglidar%beta_tot(ip,ii,k)=missing_val
    119126       endif
    120127     enddo
    … …  
    126133    do ip = 1,Npoints
    127134     if(stlidar%cldlayer(ip,k).eq.R_UNDEF)then
    128       stlidar%cldlayer(ip,k)=Cosp_fill_value
     135       stlidar%cldlayer(ip,k)=missing_val
    129136     endif
    130137    enddo
    … …  
    133140! AI 11 / 2015
    134141
    135    where(stlidar%parasolrefl == R_UNDEF) stlidar%parasolrefl = 0.0
    136    where(stlidar%lidarcldtmp == R_UNDEF) stlidar%lidarcldtmp = 0.0
    137    where(stlidar%cldlayerphase == R_UNDEF) stlidar%cldlayerphase = 0.0
    138    where(stlidar%lidarcldphase == R_UNDEF) stlidar%lidarcldphase = 0.0
    139    where(stlidar%lidarcldtmp == R_UNDEF) stlidar%lidarcldtmp = 0.0
     142   where(stlidar%parasolrefl == R_UNDEF) stlidar%parasolrefl = missing_val
     143   where(stlidar%lidarcldtmp == R_UNDEF) stlidar%lidarcldtmp = missing_val
     144   where(stlidar%cldlayerphase == R_UNDEF) stlidar%cldlayerphase = missing_val
     145   where(stlidar%lidarcldphase == R_UNDEF) stlidar%lidarcldphase = missing_val
     146   where(stlidar%lidarcldtmp == R_UNDEF) stlidar%lidarcldtmp = missing_val
    140147   
    141148
    … …  
    169176   CALL histwrite3d_cosp(o_clcalipsotmpun,stlidar%lidarcldtmp(:,:,4),nverttemp)
    170177
     178#ifdef CPP_XIOS
     179   CALL histwrite4d_cosp(o_cfad_lidarsr532,stlidar%cfad_sr)
     180#else
    171181   do icl=1,SR_BINS
    172182      CALL histwrite3d_cosp(o_cfad_lidarsr532,stlidar%cfad_sr(:,icl,:),nvert,icl)
    173183   enddo
    174 
    175    CALL histwrite3d_cosp(o_parasol_refl,stlidar%parasolrefl,nvertp)
     184     CALL histwrite3d_cosp(o_parasol_refl,stlidar%parasolrefl,nvertp)
     185#endif
    176186
    177187   do k=1,PARASOL_NREFL
    … …  
    182192        Ncref(ip,k) = 1.
    183193     else
    184         parasolcrefl(ip,k)=0.
     194        parasolcrefl(ip,k)=missing_val
    185195        Ncref(ip,k) = 0.
    186196     endif
    … …  
    190200   CALL histwrite3d_cosp(o_parasol_crefl,parasolcrefl,nvertp)
    191201
     202#ifdef CPP_XIOS
     203   CALL histwrite4d_cosp(o_atb532,sglidar%beta_tot)
     204#else
    192205   do icl=1,Ncolumns
    193206      CALL histwrite3d_cosp(o_atb532,sglidar%beta_tot(:,icl,:),nvertmcosp,icl)
    194207   enddo
     208#endif
     209
    195210   CALL histwrite3d_cosp(o_beta_mol532,sglidar%beta_mol,nvertmcosp)
    196211 endif !Lidar
    197212
    198213 if (cfg%Lradar_sim) then
     214
     215#ifdef CPP_XIOS
     216   CALL histwrite4d_cosp(o_dbze94,sgradar%Ze_tot)
     217   CALL histwrite4d_cosp(o_cfadDbze94,stradar%cfad_ze)
     218#else
    199219   do icl=1,Ncolumns
    200220      CALL histwrite3d_cosp(o_dbze94,sgradar%Ze_tot(:,icl,:),nvertmcosp,icl)
    … …  
    203223    CALL histwrite3d_cosp(o_cfadDbze94,stradar%cfad_ze(:,icl,:),nvert,icl)
    204224   enddo
     225#endif
    205226 endif
    206227
    207228 if (cfg%Llidar_sim .and. cfg%Lradar_sim) then
    208229   where(stradar%lidar_only_freq_cloud == R_UNDEF) &
    209                            stradar%lidar_only_freq_cloud = 0.0
     230                           stradar%lidar_only_freq_cloud = missing_val
    210231   CALL histwrite3d_cosp(o_clcalipso2,stradar%lidar_only_freq_cloud,nvert)
    211232   where(stradar%radar_lidar_tcc == R_UNDEF) &
    212                            stradar%radar_lidar_tcc = 0.0
     233                           stradar%radar_lidar_tcc = missing_val
    213234   CALL histwrite2d_cosp(o_cltlidarradar,stradar%radar_lidar_tcc)
    214235 endif
    … …  
    219240   do ip = 1,Npoints
    220241    if(isccp%totalcldarea(ip).eq.R_UNDEF)then
    221       isccp%totalcldarea(ip)=Cosp_fill_value
     242      isccp%totalcldarea(ip)=missing_val
    222243    endif
    223244    if(isccp%meanptop(ip).eq.R_UNDEF)then
    224       isccp%meanptop(ip)=Cosp_fill_value
     245      isccp%meanptop(ip)=missing_val
    225246    endif
    226247    if(isccp%meantaucld(ip).eq.R_UNDEF)then
    227       isccp%meantaucld(ip)=Cosp_fill_value
     248      isccp%meantaucld(ip)=missing_val
    228249    endif
    229250    if(isccp%meanalbedocld(ip).eq.R_UNDEF)then
    230       isccp%meanalbedocld(ip)=Cosp_fill_value
     251     isccp%meanalbedocld(ip)=missing_val
    231252    endif
    232253    if(isccp%meantb(ip).eq.R_UNDEF)then
    233       isccp%meantb(ip)=Cosp_fill_value
     254     isccp%meantb(ip)=missing_val
    234255    endif
    235256    if(isccp%meantbclr(ip).eq.R_UNDEF)then
    236       isccp%meantbclr(ip)=Cosp_fill_value
     257       isccp%meantbclr(ip)=missing_val
    237258    endif
    238259
    … …  
    240261     do ii=1,7
    241262     if(isccp%fq_isccp(ip,ii,k).eq.R_UNDEF)then
    242       isccp%fq_isccp(ip,ii,k)=Cosp_fill_value
     263      isccp%fq_isccp(ip,ii,k)=missing_val
    243264     endif
    244265     enddo
    … …  
    247268    do ii=1,Ncolumns
    248269     if(isccp%boxtau(ip,ii).eq.R_UNDEF)then
    249        isccp%boxtau(ip,ii)=Cosp_fill_value
     270       isccp%boxtau(ip,ii)=missing_val
    250271     endif
    251272    enddo
    … …  
    253274    do ii=1,Ncolumns
    254275     if(isccp%boxptop(ip,ii).eq.R_UNDEF)then
    255        isccp%boxptop(ip,ii)=Cosp_fill_value
     276      isccp%boxptop(ip,ii)=missing_val
    256277     endif
    257278    enddo
    … …  
    259280
    260281   CALL histwrite2d_cosp(o_sunlit,gbx%sunlit)
     282#ifdef CPP_XIOS
     283   CALL histwrite4d_cosp(o_clisccp2,isccp%fq_isccp)
     284#else
    261285   do icl=1,7
    262286   CALL histwrite3d_cosp(o_clisccp2,isccp%fq_isccp(:,icl,:),nvertisccp,icl)
    263287   enddo
     288#endif
    264289   CALL histwrite3d_cosp(o_boxtauisccp,isccp%boxtau,nvertcol)
    265290   CALL histwrite3d_cosp(o_boxptopisccp,isccp%boxptop,nvertcol)
    … …  
    278303       do k=1,MISR_N_CTH
    279304        if(misr%fq_MISR(ip,ii,k).eq.R_UNDEF)then
    280               misr%fq_MISR(ip,ii,k)=Cosp_fill_value
     305            misr%fq_MISR(ip,ii,k)=missing_val
    281306        endif
    282307       enddo
    … …  
    284309   enddo
    285310
     311#ifdef CPP_XIOS
     312   CALL histwrite4d_cosp(o_clMISR,misr%fq_MISR)
     313#else
    286314   do icl=1,7
    287315      CALL histwrite3d_cosp(o_clMISR,misr%fq_MISR(:,icl,:),nvertmisr,icl)
    288316   enddo
     317#endif
    289318 endif
    290319
    … …  
    294323  do ip=1,Npoints
    295324    if(modis%Cloud_Fraction_Low_Mean(ip).eq.R_UNDEF)then
    296        modis%Cloud_Fraction_Low_Mean(ip)=Cosp_fill_value
     325      modis%Cloud_Fraction_Low_Mean(ip)=missing_val
    297326    endif
    298327    if(modis%Cloud_Fraction_High_Mean(ip).eq.R_UNDEF)then
    299        modis%Cloud_Fraction_High_Mean(ip)=Cosp_fill_value
     328       modis%Cloud_Fraction_High_Mean(ip)=missing_val
    300329    endif
    301330    if(modis%Cloud_Fraction_Mid_Mean(ip).eq.R_UNDEF)then
    302        modis%Cloud_Fraction_Mid_Mean(ip)=Cosp_fill_value
     331       modis%Cloud_Fraction_Mid_Mean(ip)=missing_val
    303332    endif
    304333    if(modis%Cloud_Fraction_Total_Mean(ip).eq.R_UNDEF)then
    305        modis%Cloud_Fraction_Total_Mean(ip)=Cosp_fill_value
     334       modis%Cloud_Fraction_Total_Mean(ip)=missing_val
    306335    endif
    307336    if(modis%Cloud_Fraction_Water_Mean(ip).eq.R_UNDEF)then
    308        modis%Cloud_Fraction_Water_Mean(ip)=Cosp_fill_value
     337       modis%Cloud_Fraction_Water_Mean(ip)=missing_val
    309338    endif
    310339    if(modis%Cloud_Fraction_Ice_Mean(ip).eq.R_UNDEF)then
    311        modis%Cloud_Fraction_Ice_Mean(ip)=Cosp_fill_value
     340       modis%Cloud_Fraction_Ice_Mean(ip)=missing_val
    312341    endif
    313342    if(modis%Optical_Thickness_Total_Mean(ip).eq.R_UNDEF)then
    314        modis%Optical_Thickness_Total_Mean(ip)=Cosp_fill_value
     343       modis%Optical_Thickness_Total_Mean(ip)=missing_val
    315344    endif
    316345    if(modis%Optical_Thickness_Water_Mean(ip).eq.R_UNDEF)then
    317        modis%Optical_Thickness_Water_Mean(ip)=Cosp_fill_value
     346       modis%Optical_Thickness_Water_Mean(ip)=missing_val
    318347    endif
    319348    if(modis%Optical_Thickness_Ice_Mean(ip).eq.R_UNDEF)then
    320        modis%Optical_Thickness_Ice_Mean(ip)=Cosp_fill_value
     349       modis%Optical_Thickness_Ice_Mean(ip)=missing_val
    321350    endif
    322351    if(modis%Cloud_Particle_Size_Water_Mean(ip).eq.R_UNDEF)then
    323        modis%Cloud_Particle_Size_Water_Mean(ip)=Cosp_fill_value
     352       modis%Cloud_Particle_Size_Water_Mean(ip)=missing_val
    324353    endif
    325354    if(modis%Cloud_Particle_Size_Ice_Mean(ip).eq.R_UNDEF)then
    326        modis%Cloud_Particle_Size_Ice_Mean(ip)=Cosp_fill_value
     355       modis%Cloud_Particle_Size_Ice_Mean(ip)=missing_val
    327356    endif
    328357    if(modis%Cloud_Top_Pressure_Total_Mean(ip).eq.R_UNDEF)then
    329        modis%Cloud_Top_Pressure_Total_Mean(ip)=Cosp_fill_value
     358       modis%Cloud_Top_Pressure_Total_Mean(ip)=missing_val
    330359    endif
    331360    if(modis%Liquid_Water_Path_Mean(ip).eq.R_UNDEF)then
    332        modis%Liquid_Water_Path_Mean(ip)=Cosp_fill_value
     361       modis%Liquid_Water_Path_Mean(ip)=missing_val
    333362    endif
    334363    if(modis%Ice_Water_Path_Mean(ip).eq.R_UNDEF)then
    335        modis%Ice_Water_Path_Mean(ip)=Cosp_fill_value
     364       modis%Ice_Water_Path_Mean(ip)=missing_val
    336365    endif
    337366
    338367  enddo
     368
     369    where(modis%Optical_Thickness_Total_LogMean == R_UNDEF) &
     370          modis%Optical_Thickness_Total_LogMean = missing_val
     371           
     372 
     373    where(modis%Optical_Thickness_Water_LogMean == R_UNDEF) &
     374          modis%Optical_Thickness_Water_LogMean = missing_val
     375
     376    where(modis%Optical_Thickness_Ice_LogMean == R_UNDEF) &
     377          modis%Optical_Thickness_Ice_LogMean = missing_val
    339378   
    340379   CALL histwrite2d_cosp(o_cllmodis,modis%Cloud_Fraction_Low_Mean)
    … …  
    360399       do k=1,7
    361400       if(modis%Optical_Thickness_vs_Cloud_Top_Pressure(ip,ii,k).eq.R_UNDEF)then
    362           modis%Optical_Thickness_vs_Cloud_Top_Pressure(ip,ii,k)=0.
     401          modis%Optical_Thickness_vs_Cloud_Top_Pressure(ip,ii,k)=missing_val
    363402        endif
    364403       enddo
    … …  
    366405    enddo
    367406
     407#ifdef CPP_XIOS
     408   CALL histwrite4d_cosp(o_clmodis,modis%Optical_Thickness_vs_Cloud_Top_Pressure)
     409#else
    368410   do icl=1,7
    369411   CALL histwrite3d_cosp(o_clmodis, &
    370412     modis%Optical_Thickness_vs_Cloud_Top_Pressure(:,icl,:),nvertisccp,icl)           
    371413   enddo
     414#endif
    372415
    373416    where(modis%Optical_Thickness_vs_ReffIce == R_UNDEF) &
    374           modis%Optical_Thickness_vs_ReffIce = Cosp_fill_value
     417          modis%Optical_Thickness_vs_ReffIce = missing_val
    375418
    376419    where(modis%Optical_Thickness_vs_ReffLiq == R_UNDEF) &
    377           modis%Optical_Thickness_vs_ReffLiq = Cosp_fill_value
    378 
     420          modis%Optical_Thickness_vs_ReffLiq = missing_val
     421
     422#ifdef CPP_XIOS
     423!    print*,'dimension de crimodis=',size(modis%Optical_Thickness_vs_ReffIce,2),&
     424!                                    size(modis%Optical_Thickness_vs_ReffIce,3)
     425    CALL histwrite4d_cosp(o_crimodis,modis%Optical_Thickness_vs_ReffIce)
     426    CALL histwrite4d_cosp(o_crlmodis,modis%Optical_Thickness_vs_ReffLiq)
     427#else
    379428    do icl=1,7
    380429   CALL histwrite3d_cosp(o_crimodis, &
    … …  
    383432     modis%Optical_Thickness_vs_ReffLiq(:,icl,:),nvertReffLiq,icl)
    384433    enddo
     434#endif
    385435 endif
    386436
    … …  
    784834  END SUBROUTINE histwrite3d_cosp
    785835
     836! ug NOUVELLE VERSION DES WRITE AVEC LA BOUCLE DO RENTREE
     837! AI sept 2013
     838  SUBROUTINE histwrite4d_cosp(var, field)
     839  USE dimphy
     840  USE mod_phys_lmdz_para
     841  USE ioipsl
     842  use iophy
     843  USE mod_grid_phy_lmdz, ONLY: nbp_lon
     844  USE print_control_mod, ONLY: lunout,prt_level
     845
     846#ifdef CPP_XIOS
     847  USE xios, only: xios_send_field
     848#endif
     849
     850
     851  IMPLICIT NONE
     852  INCLUDE 'clesphys.h'
     853
     854    TYPE(ctrl_outcosp), INTENT(IN)    :: var
     855    REAL, DIMENSION(:,:,:), INTENT(IN)  :: field ! --> field(klon,:)
     856
     857    INTEGER :: iff, k
     858
     859    REAL,DIMENSION(klon_mpi,SIZE(field,2),SIZE(field,3)) :: buffer_omp
     860    REAL :: field4d(nbp_lon,jj_nb,SIZE(field,2),SIZE(field,3))
     861    INTEGER :: ip, n, nlev, nlev2
     862    INTEGER, ALLOCATABLE, DIMENSION(:) :: index4d
     863    CHARACTER(LEN=20) ::  nomi, nom
     864
     865  IF (prt_level >= 9) write(lunout,*)'Begin histrwrite4d ',var%name
     866
     867  IF(cosp_varsdefined) THEN
     868    !Et sinon on.... écrit
     869    IF (SIZE(field,1)/=klon) &
     870   CALL abort_physic('iophy::histwrite3d','Field first DIMENSION not equal to klon',1)           
     871
     872    nlev=SIZE(field,2)
     873    nlev2=SIZE(field,3)
     874    CALL Gather_omp(field,buffer_omp)
     875!$OMP MASTER
     876    CALL grid1Dto2D_mpi(buffer_omp,field4d)
     877
     878#ifdef CPP_XIOS
     879    IF (ok_all_xml) THEN
     880     CALL xios_send_field(var%name, Field4d(:,:,1:nlev,1:nlev2))
     881     IF (prt_level >= 1) WRITE(lunout,*)'xios_send_field ',var%name
     882    ENDIF
     883#endif
     884
     885!$OMP END MASTER   
     886  ENDIF ! vars_defined
     887  IF (prt_level >= 9) write(lunout,*)'End histrwrite4d_cosp ',nom
     888  END SUBROUTINE histwrite4d_cosp
     889
    786890  SUBROUTINE conf_cospoutputs(nam_var,cles_var)
    787891!!! Lecture des noms et cles de sortie des variables dans config.def

  • LMDZ5/branches/testing/libf/phylmd/cosp/modis_simulator.F90

    r2720 r2839  
    532532
    533533    ! ########################################################################################
    534     ! Normalize pixel counts to fraction. The first three cloud fractions have been set to -1
    535     ! in cloud-free areas, so set those places to 0.
    536     ! ########################################################################################
    537     Cloud_Fraction_High_Mean(1:nPoints)  = Cloud_Fraction_High_Mean(1:nPoints) /nSubcols
    538     Cloud_Fraction_Mid_Mean(1:nPoints)   = Cloud_Fraction_Mid_Mean(1:nPoints)  /nSubcols
    539     Cloud_Fraction_Low_Mean(1:nPoints)   = Cloud_Fraction_Low_Mean(1:nPoints)  /nSubcols
    540 
     534    ! Normalize pixel counts to fraction.
     535    ! ########################################################################################
     536    Cloud_Fraction_High_Mean(1:nPoints)  = Cloud_Fraction_High_Mean(1:nPoints)  /nSubcols
     537    Cloud_Fraction_Mid_Mean(1:nPoints)   = Cloud_Fraction_Mid_Mean(1:nPoints)   /nSubcols
     538    Cloud_Fraction_Low_Mean(1:nPoints)   = Cloud_Fraction_Low_Mean(1:nPoints)   /nSubcols
     539    Cloud_Fraction_Total_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) /nSubcols
     540    Cloud_Fraction_Ice_Mean(1:nPoints)   = Cloud_Fraction_Ice_Mean(1:nPoints)   /nSubcols
     541    Cloud_Fraction_Water_Mean(1:nPoints) = Cloud_Fraction_Water_Mean(1:nPoints) /nSubcols
     542   
    541543    ! ########################################################################################
    542544    ! Set clear-scenes to undefined

  • LMDZ5/branches/testing/libf/phylmd/cosp/phys_cosp.F90

    r2594 r2839  
    77  subroutine phys_cosp( itap,dtime,freq_cosp, &
    88                        ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
    9                         ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, &
     9                        ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
    1010                        Nptslmdz,Nlevlmdz,lon,lat, presnivs,overlaplmdz,sunlit, &
    1111                        ref_liq,ref_ice,fracTerLic,u_wind,v_wind,phis,phi,ph,p,skt,t, &
    … …  
    130130! Declaration necessaires pour les sorties IOIPSL
    131131  integer :: ii
    132   real    :: ecrit_day,ecrit_hf,ecrit_mth
     132  real    :: ecrit_day,ecrit_hf,ecrit_mth, missing_val
    133133  logical :: ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, ok_all_xml
    134134
    … …  
    192192! Allocate memory for gridbox type
    193193!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    194         print *, 'Allocating memory for gridbox type...'
     194! AI mars 2017
     195!        print *, 'Allocating memory for gridbox type...'
     196
    195197!! AI
    196198!        call construct_cosp_gridbox(dble(itap),radar_freq,surface_radar,use_mie_tables,use_gas_abs,do_ray,melt_lay,k2, &
    … …  
    215217!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    216218
    217         print *, 'Populating input structure...'
     219!        print *, 'Populating input structure...'
    218220        gbx%longitude = lon
    219221        gbx%latitude = lat
    … …  
    299301        ! Define new vertical grid
    300302!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    301         print *, 'Defining new vertical grid...'
     303!        print *, 'Defining new vertical grid...'
    302304        call construct_cosp_vgrid(gbx,Nlr,use_vgrid,csat_vgrid,vgrid)
    303305
    … …  
    305307       ! Allocate memory for other types
    306308!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    307         print *, 'Allocating memory for other types...'
     309!        print *, 'Allocating memory for other types...'
    308310        call construct_cosp_subgrid(Npoints, Ncolumns, Nlevels, sgx)
    309311        call construct_cosp_sgradar(cfg,Npoints,Ncolumns,Nlevels,N_HYDRO,sgradar)
    … …  
    325327        call cosp_output_open(Nlevlmdz, Ncolumns, presnivs, dtime, freq_cosp, &
    326328                              ok_mensuelCOSP, ok_journeCOSP, ok_hfCOSP, ok_all_xml, &
    327                               ecrit_mth, ecrit_day, ecrit_hf, use_vgrid, vgrid)
     329                              ecrit_mth, ecrit_day, ecrit_hf, use_vgrid, vgrid, stlidar)
    328330      !$OMP END MASTER
    329331      !$OMP BARRIER
    … …  
    333335        ! Call simulator
    334336!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    335         print *, 'Calling simulator...'
     337!        print *, 'Calling simulator...'
    336338!! AI
    337339!        call cosp(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,stradar,stlidar)
    … …  
    347349!!!!!!!!!!!!!!!!!! Ecreture des sorties Cosp !!!!!!!!!!!!!!r!!!!!!:!!!!!
    348350
    349        print *, 'Calling write output'
    350         call cosp_output_write(Nlevlmdz, Npoints, Ncolumns, itap, dtime, freq_COSP, &
     351!       print *, 'Calling write output'
     352        call cosp_output_write(Nlevlmdz, Npoints, Ncolumns, itap, dtime, freq_COSP, missing_val, &
    351353                               cfg, gbx, vgrid, sglidar, sgradar, stlidar, stradar, &
    352354                               isccp, misr, modis)
    … …  
    355357        ! Deallocate memory in derived types
    356358!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    357         print *, 'Deallocating memory...'
     359!        print *, 'Deallocating memory...'
    358360        call free_cosp_gridbox(gbx)
    359361        call free_cosp_subgrid(sgx)

  • LMDZ5/branches/testing/libf/phylmd/cv3_routines.F90

    r2787 r2839  
    281281                    wghti, tnk, thnk, qnk, qsnk, unk, vnk, &
    282282                    cpnk, hnk, nk, icb, icbmax, iflag, gznk, plcl)
     283
     284  USE mod_phys_lmdz_transfert_para, ONLY : bcast
    283285  IMPLICIT NONE
    284286
    … …  
    349351    END IF
    350352    PRINT *, ' ok_new_feed: ', ok_new_feed
    351     first = .FALSE.
    352353!$OMP END MASTER
     354    call bcast(ok_new_feed)
     355    first = .FALSE.   
    353356  END IF
    354357!jyg>

  • LMDZ5/branches/testing/libf/phylmd/cv3p1_closure.F90

    r2542 r2839  
    543543      ELSE IF (flag_wb==2) THEN
    544544        wbeff(il) = wbmax*(0.01*(ph(il,1)-plfc(il)))**2
     545      ELSE ! Option provisoire ou le iflag_wb/10 est considere comme une vitesse
     546        wbeff(il) = flag_wb*0.01+wbmax/(1.+500./(ph(il,1)-plfc(il)))
    545547      END IF
    546548    END IF

  • LMDZ5/branches/testing/libf/phylmd/dyn1d/1DUTILS.h

    r2720 r2839  
    27832783         hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
    27842784         vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
     2785         dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
    27852786     
    27862787         else !play>plev_prof_cas(1)
    … …  
    28092810         hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
    28102811         vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
     2812         dtrad_mod_cas(l)= frac1*dtrad_prof_cas(k1) - frac2*dtrad_prof_cas(k2)
    28112813
    28122814         endif ! play.le.plev_prof_cas(1)
    … …  
    28372839         hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                            !jyg
    28382840         vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                            !jyg
     2841         dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact                      !jyg
    28392842 
    28402843        endif ! play
    … …  
    51625165         hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
    51635166         vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
     5167         dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
    51645168     
    51655169         else !play>plev_prof_cas(1)
    … …  
    51985202         hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
    51995203         vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
     5204         dtrad_mod_cas(l)= frac1*dtrad_prof_cas(k1) - frac2*dtrad_prof_cas(k2)
    52005205
    52015206         endif ! play.le.plev_prof_cas(1)
    … …  
    52345239         hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                     !jyg
    52355240         vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                     !jyg
     5241         dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact               !jyg
    52365242 
    52375243        endif ! play

  • LMDZ5/branches/testing/libf/phylmd/dyn1d/lmdz1d.F90

    r2720 r2839  
    4040   USE physiq_mod, ONLY: physiq
    4141   USE comvert_mod, ONLY: presnivs, ap, bp, dpres,nivsig, nivsigs, pa, &
    42                           preff
     42                          preff, aps, bps, pseudoalt, scaleheight
    4343   USE temps_mod, ONLY: annee_ref, calend, day_end, day_ini, day_ref, &
    4444                        itau_dyn, itau_phy, start_time
    … …  
    634634        call disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
    635635        print *,'On utilise disvert0'
     636        aps(1:llm)=0.5*(ap(1:llm)+ap(2:llm+1))
     637        bps(1:llm)=0.5*(bp(1:llm)+bp(2:llm+1))
     638        scaleheight=8.
     639        pseudoalt(1:llm)=-scaleheight*log(presnivs(1:llm)/preff)
    636640      ELSE
    637641        call disvert()
    … …  
    640644!       Dans ce cas, on lit ap,bp dans le fichier hybrid.txt
    641645      ENDIF
    642       ! initialize ap,bp, etc. in vertical_layers_mod
     646
    643647      sig_s=presnivs/preff
    644648      plev =ap+bp*psurf

  • LMDZ5/branches/testing/libf/phylmd/ener_conserv.F90

    r2408 r2839  
    153153
    154154      bils_ec(:)=0.
     155      bils_ech(:)=0.
    155156      bils_tke(:)=0.
    156157      bils_diss(:)=0.

  • LMDZ5/branches/testing/libf/phylmd/fisrtilp.F90

    r2720 r2839  
    1717  USE cloudth_mod
    1818  USE ioipsl_getin_p_mod, ONLY : getin_p
     19  USE phys_local_var_mod, ONLY: ql_seri,qs_seri
     20  ! flag to include modifications to ensure energy conservation (if flag >0)
     21  USE add_phys_tend_mod, only : fl_cor_ebil
    1922  IMPLICIT none
    2023  !======================================================================
    … …  
    2629  !             et ilp (il pleut, L. Li)
    2730  ! Principales parties:
     31  ! P0> Thermalisation des precipitations venant de la couche du dessus
    2832  ! P1> Evaporation de la precipitation (qui vient du niveau k+1)
    2933  ! P2> Formation du nuage (en k)
     34  ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
     35  ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
     36  ! les valeurs de T et Q initiales
     37  ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
     38  ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
     39  ! P2.B> Nuage "tout ou rien"
     40  ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
    3041  ! P3> Formation de la precipitation (en k)
    3142  !======================================================================
     43  ! JLD:
     44  ! * Routine probablement fausse (au moins incoherente) si thermcep = .false.
     45  ! * fl_cor_ebil doit etre > 0 ;
     46  !   fl_cor_ebil= 0 pour reproduire anciens bugs
    3247  !======================================================================
    3348  include "YOMCST.h"
    … …  
    3853  ! Principaux inputs:
    3954  !
    40   REAL dtime ! intervalle du temps (s)
    41   REAL paprs(klon,klev+1) ! pression a inter-couche
    42   REAL pplay(klon,klev) ! pression au milieu de couche
    43   REAL t(klon,klev) ! temperature (K)
    44   REAL q(klon,klev) ! humidite specifique (kg/kg)
     55  REAL, INTENT(IN)                              :: dtime  ! intervalle du temps (s)
     56  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs  ! pression a inter-couche
     57  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay  ! pression au milieu de couche
     58  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: t      ! temperature (K)
     59  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: q      ! humidite specifique (kg/kg)
     60  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv ! points ou le schema de conv. prof. est actif
     61  INTEGER,                         INTENT(IN)   :: iflag_cld_th
     62  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo
     63  !
     64  ! Inputs lies aux thermiques
     65  !
     66  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ztv
     67  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zqta, fraca
     68  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zpspsk, ztla
     69  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zthl
     70  !
     71  !  Input/output
     72  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratqs  ! determine la largeur de distribution de vapeur
    4573  !
    4674  ! Principaux outputs:
    4775  !
    48   REAL d_t(klon,klev) ! incrementation de la temperature (K)
    49   REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
    50   REAL d_ql(klon,klev) ! incrementation de l'eau liquide
    51   REAL d_qi(klon,klev) ! incrementation de l'eau glace
    52   REAL rneb(klon,klev) ! fraction nuageuse
    53   REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
    54   REAL rhcl(klon,klev) ! humidite relative en ciel clair
    55   REAL rain(klon) ! pluies (mm/s)
    56   REAL snow(klon) ! neige (mm/s)
    57   REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
    58   REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
    59   !
    60   ! Autres arguments
    61   !
    62   REAL ztv(klon,klev)
    63   REAL zqta(klon,klev),fraca(klon,klev)
    64   REAL sigma1(klon,klev),sigma2(klon,klev)
    65   REAL qltot(klon,klev),ctot(klon,klev)
    66   REAL zpspsk(klon,klev),ztla(klon,klev)
    67   REAL zthl(klon,klev)
    68   REAL ztfondue, qsl, qsi
    69 
    70   logical lognormale(klon)
    71   logical ice_thermo
     76  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t  ! incrementation de la temperature (K)
     77  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q  ! incrementation de la vapeur d'eau
     78  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql ! incrementation de l'eau liquide
     79  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi ! incrementation de l'eau glace
     80  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb ! fraction nuageuse
     81  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radliq ! eau liquide utilisee dans rayonnements
     82  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl ! humidite relative en ciel clair
     83  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain
     84  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow
     85  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl
     86  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl
    7287
    7388  !AA
    7489  ! Coeffients de fraction lessivee : pour OFF-LINE
    7590  !
    76   REAL pfrac_nucl(klon,klev)
    77   REAL pfrac_1nucl(klon,klev)
    78   REAL pfrac_impa(klon,klev)
     91  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfrac_nucl
     92  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfrac_1nucl
     93  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfrac_impa
    7994  !
    8095  ! Fraction d'aerosols lessivee par impaction et par nucleation
    8196  ! POur ON-LINE
    8297  !
    83   REAL frac_impa(klon,klev)
    84   REAL frac_nucl(klon,klev)
    85   real zct      ,zcl
     98  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa
     99  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl
    86100  !AA
     101  ! --------------------------------------------------------------------------------
    87102  !
    88103  ! Options du programme:
    … …  
    92107
    93108  INTEGER ninter ! sous-intervals pour la precipitation
    94   INTEGER ncoreczq 
    95   INTEGER iflag_cld_th
    96   INTEGER iflag_ice_thermo
    97109  PARAMETER (ninter=5)
    98110  LOGICAL evap_prec ! evaporation de la pluie
    99111  PARAMETER (evap_prec=.TRUE.)
    100   REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
    101   logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
     112  !
     113  LOGICAL cpartiel ! condensation partielle
     114  PARAMETER (cpartiel=.TRUE.)
     115  REAL t_coup
     116  PARAMETER (t_coup=234.0)
     117  REAL DDT0
     118  PARAMETER (DDT0=.01)
     119  REAL ztfondue
     120  PARAMETER (ztfondue=278.15)
     121  ! --------------------------------------------------------------------------------
     122  !
     123  ! Variables locales:
     124  !
     125  INTEGER i, k, n, kk
     126  REAL qsl, qsi
     127  real zct      ,zcl
     128  INTEGER ncoreczq 
     129  REAL ctot(klon,klev)
     130  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 
     131  REAL zdqsdT_raw(klon)
     132  REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
     133
     134  logical lognormale(klon)
     135  logical ice_thermo
     136  LOGICAL convergence(klon)
     137  INTEGER n_i(klon), iter
     138  REAL cste
    102139
    103140  real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
    … …  
    105142  real erf   
    106143  REAL qcloud(klon)
    107   !
    108   LOGICAL cpartiel ! condensation partielle
    109   PARAMETER (cpartiel=.TRUE.)
    110   REAL t_coup
    111   PARAMETER (t_coup=234.0)
    112   !
    113   ! Variables locales:
    114   !
    115   INTEGER i, k, n, kk
    116   REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 
    117   REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
    118   LOGICAL convergence(klon)
    119   REAL DDT0
    120   PARAMETER (DDT0=.01)
    121   INTEGER n_i(klon), iter
    122   REAL cste
    123144 
    124145  REAL zrfl(klon), zrfln(klon), zqev, zqevt
    … …  
    134155  REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
    135156  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
    136   REAL zmelt, zpluie, zice, zcondold
    137   PARAMETER (ztfondue=278.15)
     157  REAL zmelt, zpluie, zice
    138158  REAL dzfice(klon)
    139159  REAL zsolid
    140160!!!!
    141161!  Variables pour Bergeron
    142   REAL zcp, coef1, DeltaT
     162  REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
    143163  REAL zqpreci(klon), zqprecl(klon)
     164! Variable pour conservation enegie des precipitations
     165  REAL zmqc(klon)
    144166  !
    145167  LOGICAL appel1er
    … …  
    170192  REAL   beta(klon,klev) ! taux de conversion de l'eau cond
    171193! RomP <<<
    172   REAL zmair, zcpair, zcpeau
     194  REAL zmair(klon), zcpair, zcpeau
    173195  !     Pour la conversion eau-neige
    174196  REAL zlh_solid(klon), zm_solid
    … …  
    180202                     ! (Heymsfield & Donner, 1990)
    181203  REAL zzz
     204
    182205  include "YOETHF.h"
    183206  include "FCTTRE.h"
    … …  
    312335     ENDDO
    313336     !
     337     ! ----------------------------------------------------------------
     338     ! P0> Thermalisation des precipitations venant de la couche du dessus
     339     ! ----------------------------------------------------------------
    314340     ! Calculer la varition de temp. de l'air du a la chaleur sensible
    315      ! transporter par la pluie.
    316      ! Il resterait a rajouter cet effet de la chaleur sensible sur les
    317      ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
    318      ! surface.
     341     ! transporter par la pluie. On thermalise la pluie avec l'air de la couche.
     342     ! Cette quantite de pluie qui est thermalisee, et devra continue a l'etre lors
     343     ! des differentes transformations thermodynamiques. Cette masse d'eau doit
     344     ! donc etre ajoute a l'humidite de la couche lorsque l'on calcule la variation
     345     ! de l'enthalpie  de la couche avec la temperature
     346     ! Variables calculees ou modifiees:
     347     !   -  zt: temperature de la cocuhe
     348     !   - zmqc: masse de precip qui doit etre thermalisee
    319349     !
    320350     IF(k.LE.klevm1) THEN         
    321351        DO i = 1, klon
    322352           !IM
    323            zmair=(paprs(i,k)-paprs(i,k+1))/RG
     353           zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
     354           ! il n'y a pas encore d'eau liquide ni glace dans la maiille, donc zq suffit
    324355           zcpair=RCPD*(1.0+RVTMP2*zq(i))
    325356           zcpeau=RCPD*RVTMP2
     357         if (fl_cor_ebil .GT. 0) then
     358           ! zmqc: masse de precip qui doit etre thermalisee avec l'air de la couche atm
     359           ! pour s'assurer que la precip arrivant au sol aura bien la temperature de la
     360           ! derniere couche
     361           zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
     362           ! t(i,k+1)+d_t(i,k+1): nouvelle temp de la couche au dessus
     363           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
     364                 / (zcpair + zmqc(i)*zcpeau)
     365         else ! si on maintient les anciennes erreurs
    326366           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
    327                 + zmair*zcpair*zt(i) ) &
    328                 / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
    329            !     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
     367                + zmair(i)*zcpair*zt(i) ) &
     368                / (zmair(i)*zcpair + zrfl(i)*dtime*zcpeau)
     369         end if
    330370        ENDDO
    331      ENDIF
    332      ! ----------------------------------------------------------------
    333      ! P1> Debut evaporation de la precipitation
    334      ! ----------------------------------------------------------------
     371     ENDIF ! end IF(k.LE.klevm1)
     372     !
     373     ! ----------------------------------------------------------------
     374     ! P1> Calcul de l'evaporation de la precipitation
     375     ! ----------------------------------------------------------------
     376     ! On evapore une partie des precipitations venant de la maille du dessus.
     377     ! On calcule l'evaporation et la sublimation des precipitations, jusqu'a
     378     ! ce que la fraction de cette couche qui est sous le nuage soit saturee.
     379     ! Variables calculees ou modifiees:
     380     !   - zrfl et zifl: flux de precip liquide et glace
     381     !   - zq, zt: humidite et temperature de la cocuhe
     382     !   - zmqc: masse de precip qui doit etre thermalisee
     383     !
    335384     IF (evap_prec) THEN
    336385        DO i = 1, klon
    337 !AJ<
    338 !!           IF (zrfl(i) .GT.0.) THEN
     386!          S'il y a des precipitations
    339387           IF (zrfl(i)+zifl(i).GT.0.) THEN
    340 !>AJ
    341388              ! Calcul du qsat
    342389              IF (thermcep) THEN
    … …  
    356403        ENDDO
    357404!AJ<
     405
    358406       IF (.NOT. ice_thermo) THEN
    359407        DO i = 1, klon
    360 !AJ<
    361 !!           IF (zrfl(i) .GT.0.) THEN
     408!          S'il y a des precipitations
    362409           IF (zrfl(i)+zifl(i).GT.0.) THEN
    363 !>AJ
    364410                ! Evap max pour ne pas saturer la fraction sous le nuage
     411                ! Evap max jusqu'à atteindre la saturation dans la partie
     412                ! de la maille qui est sous le nuage de la couche du dessus
     413                !!! On ne tient compte de cette fraction que sous une seule
     414                !!! couche sous le nuage
    365415                zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
     416             ! Ajout de la prise en compte des precip a thermiser
     417             ! avec petite reecriture
     418             if  (fl_cor_ebil .GT. 0) then ! nouveau
     419                ! Calcul de l'evaporation du flux de precip herite
     420                !   d'au-dessus
     421                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
     422                     * zmair(i)/pplay(i,k)*zt(i)*RD
     423                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * dtime/zmair(i)
     424
     425                ! Seuil pour ne pas saturer la fraction sous le nuage
     426                zqev = MIN (zqev, zqevt)
     427                ! Nouveau flux de precip
     428                zrfln(i) = zrfl(i) - zqev*zmair(i)/dtime
     429                ! Aucun flux liquide pour T < t_coup, on reevapore tout.
     430                IF (zt(i) .LT. t_coup.and.reevap_ice) THEN
     431                  zrfln(i)=0.
     432                  zqev = (zrfl(i)-zrfln(i))/zmair(i)*dtime
     433                END IF
     434                ! Nouvelle vapeur
     435                zq(i) = zq(i) + zqev
     436                zmqc(i) = zmqc(i)-zqev
     437                ! Nouvelle temperature (chaleur latente)
     438                zt(i) = zt(i) - zqev &
     439                     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
     440!!JLD debut de partie a supprimer a terme
     441            else ! if  (fl_cor_ebil .GT. 0)
    366442                ! Calcul de l'evaporation du flux de precip herite
    367443                !   d'au-dessus
    … …  
    384460                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
    385461                     * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
     462              end if ! if  (fl_cor_ebil .GT. 0)
     463!!JLD fin de partie a supprimer a terme
    386464                zrfl(i) = zrfln(i)
    387465                zifl(i) = 0.
    … …  
    390468!
    391469       ELSE ! (.NOT. ice_thermo)
    392 !
     470!      ================================
     471!      Avec thermodynamique de la glace
     472!      ================================
    393473        DO i = 1, klon
    394474!AJ<
    395 !!           IF (zrfl(i) .GT.0.) THEN
    396            IF (zrfl(i)+zifl(i).GT.0.) THEN
    397 !>AJ
    398      !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    399      ! Modification de l'evaporation avec la glace
    400      ! Differentiation entre precipitation liquide et solide
    401      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     475!        S'il y a des precipitations
     476         IF (zrfl(i)+zifl(i).GT.0.) THEN
    402477     
    403      ! Evap max pour ne pas saturer la fraction sous le nuage
     478        ! Evap max pour ne pas saturer la fraction sous le nuage
    404479         zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
    405      !    zqev0 = MAX (0.0, zqs(i)-zq(i) )
    406 
    407      !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    408      ! On differencie qsat pour l'eau et la glace
    409      ! Si zdelta=1. --> glace
    410      ! Si zdelta=0. --> eau liquide
    411      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     480
     481         !JAM
     482         ! On differencie qsat pour l'eau et la glace
     483         ! Si zdelta=1. --> glace
     484         ! Si zdelta=0. --> eau liquide
    412485       
    413486         ! Calcul du qsat par rapport a l'eau liquide
    … …  
    417490         qsl= qsl*zcor
    418491         
    419          ! Calcul de l'evaporation du flux de precip herite
    420          !   d'au-dessus
     492         ! Calcul de l'evaporation du flux de precip venant du dessus
    421493         ! Formulation en racine du flux de precip
    422494         ! (Klemp & Wilhelmson, 1978; Sundqvist, 1988)
    … …  
    425497         zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
    426498              *RG*dtime/(paprs(i,k)-paprs(i,k+1))
    427 
    428499         
    429500         ! Calcul du qsat par rapport a la glace
    … …  
    440511              *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
    441512
    442              
    443      !JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    444      ! Verification sur l'evaporation
    445      ! On s'assure qu'on ne sature pas
    446      ! la fraction sous le nuage sinon on
    447      ! repartit zqev0 en gardant la proportion
    448      ! liquide / glace
    449      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     513        !JAM
     514        ! Limitation de l'evaporation. On s'assure qu'on ne sature pas
     515        ! la fraction de la couche sous le nuage sinon on repartit zqev0
     516        ! en conservant la proportion liquide / glace
    450517     
    451518         IF (zqevt+zqevti.GT.zqev0) THEN
    452                   zqev=zqev0*zqevt/(zqevt+zqevti)
    453                   zqevi=zqev0*zqevti/(zqevt+zqevti)
    454              
     519            zqev=zqev0*zqevt/(zqevt+zqevti)
     520            zqevi=zqev0*zqevti/(zqevt+zqevti)
    455521         ELSE
     522!JLD je ne comprends pas les lignes ci-dessous. On repartit les precips
     523!       liquides et solides meme si on ne sature pas la couche.
     524!       A mon avis, le test est inutile, et il faudrait tout remplacer par:
     525!            zqev=zqevt
     526!            zqevi=zqevti
    456527             IF (zqevt+zqevti.GT.0.) THEN
    457                   zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
    458                   zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
     528                zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
     529                zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
    459530             ELSE
    460531             zqev=0.
    … …  
    471542         zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
    472543                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
     544       if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
     545         zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
     546                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
     547         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
     548                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
     549                  * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
     550                  + (zifln(i)-zifl(i)) &
     551                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
     552                  * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
     553       else ! sans correction thermalisation des precips
    473554         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
    474555                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
    … …  
    477558                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
    478559                  * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
    479      
     560       end if
     561        ! Nouvelles vaeleurs des precips liquides et solides
    480562         zrfl(i) = zrfln(i)
    481563         zifl(i) = zifln(i)
    … …  
    486568!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
    487569           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))             ! jyg
     570           ! fraction de la precip solide qui est fondue
    488571           zmelt = MIN(MAX(zmelt,0.),1.)
    489572           ! Fusion de la glace
    490573           zrfl(i)=zrfl(i)+zmelt*zifl(i)
    491            zifl(i)=zifl(i)*(1.-zmelt)
    492 !           print*,zt(i),'octavio1'
     574           if (fl_cor_ebil .LE. 0) then
     575             ! the following line should not be here. Indeed, if zifl is modified
     576             ! now, zifl(i)*zmelt is no more the amount of ice that has melt
     577             ! and therefore the change in temperature computed below is wrong
     578             zifl(i)=zifl(i)*(1.-zmelt)
     579           end if
    493580           ! Chaleur latente de fusion
     581        if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
     582           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
     583                      *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
     584        else ! sans correction thermalisation des precips
    494585           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
    495586                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
    496 !           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
    497 !fin CR
    498 
    499 
     587        end if
     588           if (fl_cor_ebil .GT. 0) then ! correction bug, deplacement ligne precedente
     589             zifl(i)=zifl(i)*(1.-zmelt)
     590           end if
    500591
    501592           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
    … …  
    515606           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
    516607           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
     608       if  (fl_cor_ebil .GT. 0) then ! nouveau
     609           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
     610       else   
    517611           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
     612       endif
    518613           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
    519614           zqs(i) = MIN(0.5,zqs(i))
    … …  
    521616           zqs(i) = zqs(i)*zcor
    522617           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
     618           zdqsdT_raw(i) = zdqs(i)*  &
     619         &         RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
    523620        ENDDO
    524621     ELSE
    … …  
    547644     ! P2> Formation du nuage
    548645     ! ----------------------------------------------------------------
     646     ! Variables calculees:
     647     !   rneb  : fraction nuageuse
     648     !   zcond : eau condensee moyenne dans la maille.
     649     !   rhcl: humidite relative ciel-clair
     650     !   zt : temperature de la maille
     651     ! ----------------------------------------------------------------
     652     !
    549653     IF (cpartiel) THEN
    550 
    551         !        print*,'Dans partiel k=',k
     654        ! -------------------------
     655        ! P2.A> Nuage fractionnaire
     656        ! -------------------------
    552657        !
    553658        !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
    … …  
    561666        !   Version avec les raqts
    562667
     668        ! ----------------------------------------------------------------
     669        ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
     670        ! ----------------------------------------------------------------
    563671        if (iflag_pdf.eq.0) then
    564672
    … …  
    570678           enddo
    571679
    572         else
    573            !
    574            !   Version avec les nouvelles PDFs.
     680        else !  if (iflag_pdf.eq.0)
     681           ! ----------------------------------------------------------------
     682           ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
     683           ! les valeurs de T et Q initiales
     684           ! ----------------------------------------------------------------
    575685           do i=1,klon
    576686              if(zq(i).lt.1.e-15) then
    … …  
    620730           enddo
    621731
     732        ! ----------------------------------------------------------------
     733        ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
     734        ! Calcul des grandeurs nuageuses en tenant compte de l'effet de
     735        ! la condensation sur T, et donc sur qsat et sur les grandeurs nuageuses
     736        ! qui en dependent. Ce changement de temperature est provisoire, et
     737        ! la valeur definitive sera calcule plus tard.
     738        ! Variables calculees:
     739        !   rneb : nebulosite
     740        !   zcond: eau condensee en moyenne dans la maille
     741        ! note JLD: si on n'a pas de pdf lognormale, ce qui se passe ne me semble
     742        ! pas clair, il n'y a probablement pas de prise en compte de l'effet de
     743        ! T sur qsat
     744        ! ----------------------------------------------------------------
    622745
    623746!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
    … …  
    629752               do i=1,klon
    630753!                 do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
     754                 ! !! convergence = .true. tant que l'on n'a pas converge !!
     755                 ! ------------------------------
    631756                 convergence(i)=abs(DT(i)).gt.DDT0
    632757                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
    633                  Tbef(i)=Tbef(i)+DT(i)
     758                 ! si on n'a pas converge
     759                 !
     760                 ! P2.A.2.1> Calcul de la fraction nuageuse et de la quantite d'eau condensee
     761                 ! ---------------------------------------------------------------
     762                 ! Variables calculees:
     763                 ! rneb : nebulosite
     764                 ! zqn : eau condensee, dans le nuage (in cloud water content)
     765                 ! zcond: eau condensee en moyenne dans la maille
     766                 ! rhcl: humidite relative ciel-clair
     767                 !
     768                 Tbef(i)=Tbef(i)+DT(i) ! nouvelle temperature
    634769                 if (.not.ice_thermo) then   
    635770                    zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
    … …  
    641776                          zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i)))
    642777                       else
    643 !avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
     778                          !avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
    644779                          zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
    645780                       endif
    646781                    endif
    647782                 endif
    648                  ! Calcul des PDF lognormales
     783                 ! Calcul de rneb, qzn et zcond pour les PDF lognormales
    649784                 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
     785               if (fl_cor_ebil .GT. 0) then
     786                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
     787               else
    650788                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
     789               end if
    651790                 zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
    652791                 zqs(i) = MIN(0.5,zqs(i))
    … …  
    677816               enddo ! boucle en i
    678817
     818                 ! P2.A.2.2> Calcul APPROCHE de la variation de temperature DT
     819                 !         due a la condensation.
     820                 ! ---------------------------------------------------------------
     821                 ! Variables calculees:
     822                 ! DT : variation de temperature due a la condensation
     823
    679824                 if (.not. ice_thermo) then
    680 
     825                 ! --------------------------
    681826                 do i=1,klon
    682827                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
    683828
    684829                 qlbef(i)=max(0.,zqn(i)-zqs(i))
     830               if (fl_cor_ebil .GT. 0) then
     831                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
     832               else
    685833                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
     834               end if
    686835                 denom=1.+rneb(i,k)*zdqs(i)
    687836                 DT(i)=num/denom
    … …  
    690839                 enddo
    691840
    692                  else
    693                  ! Iteration pour convergence avec qsat(T)
     841                 else ! if (.not. ice_thermo)
     842                 ! --------------------------
    694843                 if (iflag_t_glace.ge.1) then
    695844                 CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
    … …  
    703852                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
    704853                    zfice(i) = zfice(i)**exposant_glace_old
    705                     dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) / (t_glace_min_old - RTT)
     854                    dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) &
     855          &                     / (t_glace_min_old - RTT)
    706856                 endif
    707857                 
    708858                 if (iflag_t_glace.ge.1) then
    709                  dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) / (t_glace_min - t_glace_max)
     859                 dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) &
     860          &                    / (t_glace_min - t_glace_max)
    710861                 endif
    711862               
    … …  
    721872
    722873                 qlbef(i)=max(0.,zqn(i)-zqs(i))
    723                  num=-Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
    724                  denom=1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
     874               if (fl_cor_ebil .GT. 0) then
     875                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
     876           &          +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
     877                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
     878                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)    &
     879           &               *qlbef(i)*dzfice(i)
     880               else
     881                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
     882           &         +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
     883                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
    725884                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))*rneb(i,k)*qlbef(i)*dzfice(i)
     885               end if
    726886                 DT(i)=num/denom
    727887                 n_i(i)=n_i(i)+1
    … …  
    732892                 endif !ice_thermo
    733893
    734 !                 endif
    735 !               enddo
    736  
    737          
    738894             enddo ! iter=1,iflag_fisrtilp_qsat+1
    739895             ! Fin d'iteration pour condensation avec variation de qsat(T)
    740896             ! -----------------------------------------------------------
    741            endif
    742 
     897           endif !  if (iflag_fisrtilp_qsat.ge.0)
     898     ! ----------------------------------------------------------------
     899     ! Fin de P2.A.2> la prise en compte du couplage entre eau condensee et T
     900     ! ----------------------------------------------------------------
    743901
    744902        endif ! iflag_pdf
    745 
    746903
    747904!        if (iflag_fisrtilp_qsat.eq.-1) then
    … …  
    771928
    772929!        ELSE
    773 
    774         ! Calcul de l'eau in-cloud (zqn),
    775         ! moyenne dans la maille (zcond),
    776         ! fraction nuageuse (rneb) et
    777         ! humidite relative ciel-clair (rhcl)
     930        ! ----------------------------------------------------------------
     931        ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
     932        ! Variables calculees:
     933        !   rneb : nebulosite
     934        !   zcond: eau condensee en moyenne dans la maille
     935        !   zq : eau vapeur dans la maille
     936        !   zt : temperature de la maille
     937        !   rhcl: humidite relative ciel-clair
     938        ! ----------------------------------------------------------------
     939        !
     940        ! Bornage de l'eau in-cloud (zqn) et de la fraction nuageuse (rneb)
     941        ! Calcule de l'eau condensee moyenne dans la maille (zcond),
     942        ! et de l'humidite relative ciel-clair (rhcl)
    778943        DO i=1,klon
    779944           IF (rneb(i,k) .LE. 0.0) THEN
    … …  
    793958        ENDDO
    794959
    795 
    796960!        ENDIF
    797961
    798962     ELSE ! de IF (cpartiel)
    799         ! Cas "tout ou rien"
     963        ! -------------------------
     964        ! P2.B> Nuage "tout ou rien"
     965        ! -------------------------
     966        ! note JLD: attention, rhcl non calcule. Ca peut avoir des consequences?
    800967        DO i = 1, klon
    801968           IF (zq(i).GT.zqs(i)) THEN
    … …  
    806973           zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
    807974        ENDDO
    808      ENDIF
    809      ! ----------------------------------------------------------------
    810      ! Fin de formation du nuage
    811      ! ----------------------------------------------------------------
     975     ENDIF ! de IF (cpartiel)
    812976     !
    813977     ! Mise a jour vapeur d'eau
     978     ! -------------------------
    814979     DO i = 1, klon
    815980        zq(i) = zq(i) - zcond(i)
    … …  
    817982     ENDDO
    818983!AJ<
    819      ! Chaleur latente apres formation nuage
     984     ! ------------------------------------
     985     ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
    820986     ! -------------------------------------
     987     ! Variable calcule:
     988     !   zt : temperature de la maille
     989     !
    821990     IF (.NOT. ice_thermo) THEN
    822991        if (iflag_fisrtilp_qsat.lt.1) then
    … …  
    826995        else if (iflag_fisrtilp_qsat.gt.0) then
    827996           DO i= 1, klon
     997    if (fl_cor_ebil .GT. 0) then
     998             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
     999    else
    8281000             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
     1001    end if
    8291002           ENDDO
    8301003        endif
    … …  
    8561029                    zfice(i) = zfice(i)**exposant_glace_old
    8571030              endif
     1031        if (fl_cor_ebil .GT. 0) then
     1032              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
     1033           &             * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
     1034                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
     1035        else
    8581036              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
    859                        +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
     1037                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
     1038        end if
    8601039           ENDDO
    8611040         endif
    … …  
    8631042       ENDIF
    8641043!>AJ
     1044
    8651045     ! ----------------------------------------------------------------
    8661046     ! P3> Formation des precipitations
    … …  
    9581138     !         zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie   , 0.0)
    9591139     !         zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice   , 0.0)
     1140!JLD : les 2 variables zoliqp et zoliqi crorresponent a des pseudo precip
     1141!      temporaires et ne doivent pas etre calcule (alors qu'elles le sont
     1142!      si iflag_bergeron <> 2
     1143!      A SUPPRIMER A TERME
    9601144              zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie  , 0.0)
    9611145              zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice  , 0.0)
    … …  
    9661150        ENDDO  ! i = 1,klon
    9671151     ENDDO     ! n = 1,ninter
     1152
    9681153     ! ----------------------------------------------------------------
    9691154     !
    … …  
    9941179             zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
    9951180             zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
     1181           if (fl_cor_ebil .GT. 0) then
     1182             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
     1183             coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
     1184!            Calcul de DT si toute les precips liquides congelent
     1185             DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
     1186!            On ne veut pas que T devienne superieur a la temp. de congelation.
     1187!            donc que Delta > RTT-zt(i
     1188             DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
     1189             zt(i)      = zt(i)      + DeltaT
     1190!            Eau vaporisee du fait de l'augmentation de T
     1191             Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
     1192!            on reajoute cette eau vaporise a la vapeur et on l'enleve des precips
     1193             zq(i) = zq(i) +  Deltaq
     1194!            Les 3 max si dessous prtotegent uniquement des erreurs d'arrondies
     1195             zcond(i)   = max( zcond(i)- Deltaq, 0. )
     1196!            precip liquide qui congele ou qui s'evapore
     1197             Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
     1198             zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
     1199!            bilan eau glacee
     1200             zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
     1201           else ! if (fl_cor_ebil .GT. 0)
     1202!            ancien calcul
    9961203             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zcond(i)))
    9971204             coef1 = RLMLT*zdqs(i)/RLVTT
    … …  
    10021209             zq(i)      = zq(i)      + zcp/RLVTT*zdqs(i)*DeltaT
    10031210             zt(i)      = zt(i)      + DeltaT
     1211           end if ! if (fl_cor_ebil .GT. 0)
    10041212           ENDIF  ! rneb(i,k) .GT. 0.0
    10051213         ENDDO
    … …  
    10211229         IF (rneb(i,k).GT.0.0) THEN
    10221230!CR on prend en compte la phase glace
    1023            if (.not.ice_thermo) then
    1024            d_ql(i,k) = zoliq(i)
    1025            d_qi(i,k) = 0.
    1026            else
     1231!JLD inutile car on ne passe jamais ici si .not.ice_thermo
     1232!           if (.not.ice_thermo) then
     1233!           d_ql(i,k) = zoliq(i)
     1234!           d_qi(i,k) = 0.
     1235!           else
    10271236           d_ql(i,k) = (1-zfice(i))*zoliq(i)
    10281237           d_qi(i,k) = zfice(i)*zoliq(i)
    1029            endif
     1238!           endif
    10301239!AJ<
    10311240           zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
    … …  
    10431252              zifl(i) = zifl(i)+zrfl(i)
    10441253              zrfl(i) = 0.
     1254           if (fl_cor_ebil .GT. 0) then
     1255              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
     1256                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
     1257           else
    10451258              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
    10461259                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))
     1260           end if
    10471261           ENDIF  ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)
    10481262!RC   
    … …  
    10941308     ! Fin de formation des precipitations
    10951309     ! ----------------------------------------------------------------
    1096      !
    10971310     !
    10981311     ! Calculer les tendances de q et de t:
    … …  
    12021415     DO i = 1, klon
    12031416        zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
    1204         zmair=(paprs(i,k)-paprs(i,k+1))/RG
     1417        zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
    12051418        zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
    1206         d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
     1419        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair(i))
    12071420     END DO
    12081421  END DO

  • LMDZ5/branches/testing/libf/phylmd/iophys.F90

    r2787 r2839  
     1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2! Interface pour ecrire en netcdf avec les routines d'enseignement
     3! iotd de Frederic Hourdin
     4!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     5
    16      subroutine iophys_ecrit(nom,lllm,titre,unite,px)
    27
    … …  
    6267      end
    6368
     69!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     70! Version avec reindexation pour appeler depuis les routines internes
     71! à la sous surface
     72!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     73
     74
     75    subroutine iophys_ecrit_index(nom,lllm,titre,unite,knon,knindex,px)
     76
     77    USE mod_phys_lmdz_para, ONLY: klon_omp
     78    USE dimphy, ONLY : klon
     79    USE mod_grid_phy_lmdz, ONLY: klon_glo
     80    IMPLICIT NONE
     81
     82! This subroutine returns the sea surface temperature already read from limit.nc
     83!
     84
     85! Arguments on input:
     86    INTEGER lllm
     87    CHARACTER (len=*) :: nom,titre,unite
     88    REAL px(klon_omp,lllm)
     89    INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
     90    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex  ! grid point number for compressed grid
     91    REAL, DIMENSION(klon,lllm) :: varout
     92
     93    INTEGER :: i,l
     94
     95    IF (klon/=klon_omp) THEN
     96      print*,'klon, klon_omp',klon,klon_omp
     97      STOP'probleme de dimension parallele'
     98    ENDIF
     99
     100    varout(1:klon,1:lllm)=0.
     101    DO l = 1, lllm
     102    DO i = 1, knon
     103       varout(knindex(i),l) = px(i,l)
     104    END DO
     105    END DO
     106    CALL iophys_ecrit(nom,lllm,titre,unite,varout)
     107
     108  END SUBROUTINE iophys_ecrit_index
     109
    64110!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    65111      SUBROUTINE iophys_ini
    … …  
    68114      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
    69115      USE dimphy, ONLY: klev
     116      USE mod_grid_phy_lmdz, ONLY: klon_glo
    70117
    71118      IMPLICIT NONE
    … …  
    88135
    89136real pi
     137INTEGER nlat_eff
    90138
    91139!   Arguments:
    … …  
    94142!$OMP MASTER
    95143    IF (is_mpi_root) THEN       
     144
     145! Bidouille pour gerer le fait que lat_reg contient deux latitudes
     146! en version uni-dimensionnelle (chose qui pourrait être résolue
     147! par ailleurs)
     148IF (klon_glo==1) THEN
     149   nlat_eff=1
     150ELSE
     151   nlat_eff=size(lat_reg)
     152ENDIF
    96153pi=2.*asin(1.)
    97154call iotd_ini('phys.nc   ', &
    98 !iim,jjp1,llm,rlonv(1:iim)*180./pi,rlatu*180./pi,presnivs)
    99 size(lon_reg),size(lat_reg),klev,lon_reg(:)*180./pi,lat_reg*180./pi,presnivs)
     155size(lon_reg),nlat_eff,klev,lon_reg(:)*180./pi,lat_reg*180./pi,presnivs)
    100156    ENDIF
    101157!$OMP END MASTER

  • LMDZ5/branches/testing/libf/phylmd/limit_read_mod.F90

    r2787 r2839  
    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/branches/testing/libf/phylmd/open_climoz_m.F90

    r1910 r2839  
    11! $Id$
    2 module open_climoz_m
     2MODULE open_climoz_m
    33
    4   implicit none
     4  USE print_control_mod, ONLY: lunout
     5  IMPLICIT NONE
    56
    6 contains
     7CONTAINS
    78
    8   subroutine open_climoz(ncid, press_in_edg)
     9!-------------------------------------------------------------------------------
     10!
     11SUBROUTINE open_climoz(ncid, press_in_cen, press_in_edg, time_in, daily, adjust)
     12!
     13!-------------------------------------------------------------------------------
     14  USE netcdf95, ONLY: nf95_open, nf95_close, nf95_gw_var, nf95_inq_varid
     15  USE netcdf,   ONLY: nf90_nowrite
     16  USE mod_phys_lmdz_mpi_data,      ONLY: is_mpi_root
     17  USE mod_phys_lmdz_mpi_transfert, ONLY: bcast_mpi
     18  USE phys_cal_mod,                ONLY: calend, year_len, year_cur
     19!-------------------------------------------------------------------------------
     20! Purpose: This procedure should be called once per "gcm" run, by a single
     21!          thread of each MPI process.
     22!          The root MPI process opens "climoz_LMDZ.nc", reads the pressure
     23! levels and the times and broadcasts them to the other processes.
     24!          We assume that, in "climoz_LMDZ.nc", the pressure levels are in hPa
     25! and in strictly ascending order.
     26!-------------------------------------------------------------------------------
     27! Arguments (OpenMP shared):
     28  INTEGER, INTENT(OUT):: ncid      !--- "climoz_LMDZ.nc" identifier
     29  REAL, POINTER :: press_in_cen(:) !--- at cells centers
     30  REAL, POINTER :: press_in_edg(:) !--- at the interfaces (pressure intervals)
     31  REAL, POINTER :: time_in(:)      !--- records times, in days since Jan. 1st
     32  LOGICAL, INTENT(IN) :: daily     !--- daily files (calendar dependent days nb)
     33  LOGICAL, INTENT(IN) :: adjust    !--- tropopause adjustement required
     34!   pressure levels press_in_cen/edg are in Pa a,d strictly ascending order.
     35!   time_in is only used for monthly files (14 records)
     36!-------------------------------------------------------------------------------
     37! Local variables:
     38  INTEGER :: varid                 !--- NetCDF variables identifier
     39  INTEGER :: nlev, ntim            !--- pressure levels and time records numbers
     40  CHARACTER(LEN=80)  :: sub
     41  CHARACTER(LEN=320) :: msg
     42!-------------------------------------------------------------------------------
     43  sub="open_climoz"
     44  WRITE(lunout,*)"Entering routine "//TRIM(sub)
    945
    10     ! This procedure should be called once per "gcm" run, by a single
    11     ! thread of each MPI process.
    12     ! The root MPI process opens "climoz_LMDZ.nc", reads the pressure
    13     ! levels and broadcasts them to the other processes.
     46  IF(is_mpi_root) THEN
    1447
    15     ! We assume that, in "climoz_LMDZ.nc", the pressure levels are in hPa
    16     ! and in strictly ascending order.
     48    !--- OPEN FILE, READ PRESSURE LEVELS AND TIME VECTOR
     49    CALL nf95_open("climoz_LMDZ.nc", nf90_nowrite, ncid)
     50    CALL nf95_inq_varid(ncid, "plev", varid)
     51    CALL nf95_gw_var(ncid, varid, press_in_cen)
     52    ! Convert from hPa to Pa because "paprs" and "pplay" are in Pa:
     53    press_in_cen = press_in_cen * 100.
     54    nlev = SIZE(press_in_cen)
     55    CALL NF95_INQ_VARID(ncID, "time", varID)
     56    CALL NF95_GW_VAR(ncid, varid, time_in)
     57    ntim = SIZE(time_in)
    1758
    18     use netcdf95, only: nf95_open, nf95_close, nf95_gw_var, nf95_inq_varid
    19     use netcdf, only: nf90_nowrite
     59    !--- BUILD EDGES OF PRESSURE INTERVALS: HALFWAY IN LOGARITHMS
     60    ALLOCATE(press_in_edg(nlev+1))
     61    press_in_edg=[0.,SQRT(press_in_cen(1:nlev-1)*press_in_cen(2:nlev)),HUGE(0.)]
    2062
    21     use mod_phys_lmdz_mpi_data, only: is_mpi_root
    22     use mod_phys_lmdz_mpi_transfert, only: bcast_mpi ! broadcast
     63    !--- CHECK RECORDS NUMBER AND DISPLAY CORRESPONDING INFORMATION
     64    IF(daily.AND.ntim/=year_len) THEN
     65      WRITE(msg,'(a,3(i4,a))')TRIM(sub)//': Expecting a daily ozone file with',&
     66     &year_len,' records (year ',year_cur,') ; found ',ntim,' instead'
     67      CALL abort_physic(sub, msg, 1)
     68    ELSE IF(ALL([360,14]/=ntim)) THEN
     69      WRITE(msg,'(a,i4,a)')TRIM(sub)//': Expecting an ozone file with 14 (mont'&
     70     &//'hly case) or 360 (old style files) records ; found ',ntim,' instead'
     71      CALL abort_physic(sub, msg, 1)
     72    ELSE
     73      IF(daily) THEN
     74         WRITE(msg,'(a,2(i4,a))')'daily file (',ntim,' days in ',year_cur,')'
     75       ELSE IF(ntim==14) THEN
     76         msg='14 records monthly file'
     77       ELSE
     78         msg='360 records files (old convention)'
     79       END IF
     80      WRITE(lunout,*)TRIM(sub)//': Using a '//TRIM(msg)
     81    END IF
    2382
    24     integer, intent(out):: ncid ! of "climoz_LMDZ.nc", OpenMP shared
     83    !--- MESSAGE ABOUT OPTIONAL STRETCHING FOR TROPOPAUSES MATCHING
     84    IF(adjust) THEN
     85      WRITE(lunout,*)TRIM(sub)//': Adjusting O3 field to match gcm tropopause.'
     86    ELSE
     87      WRITE(lunout,*)TRIM(sub)//': Interpolating O3 field directly on gcm levels.'
     88    END IF
    2589
    26     real, pointer:: press_in_edg(:)
    27     ! edges of pressure intervals for ozone climatology, in Pa, in strictly
    28     ! ascending order, OpenMP shared
     90  END IF
     91  CALL bcast_mpi(nlev)
     92  IF(.NOT.is_mpi_root) ALLOCATE(press_in_cen(nlev  )); CALL bcast_mpi(press_in_cen)
     93  IF(.NOT.is_mpi_root) ALLOCATE(press_in_edg(nlev+1)); CALL bcast_mpi(press_in_edg)
     94  CALL bcast_mpi(ntim)
     95  IF(.NOT.is_mpi_root) ALLOCATE(time_in(ntim));        CALL bcast_mpi(time_in)
    2996
    30     ! Variables local to the procedure:
     97END SUBROUTINE open_climoz
     98!
     99!-------------------------------------------------------------------------------
    31100
    32     real, pointer:: plev(:)
    33     ! (pressure levels for ozone climatology, converted to Pa, in strictly
    34     ! ascending order)
    35 
    36     integer varid ! for NetCDF
    37     integer n_plev ! number of pressure levels in the input data
    38     integer k
    39 
    40     !---------------------------------------
    41 
    42     print *, "Call sequence information: open_climoz"
    43 
    44     if (is_mpi_root) then
    45        call nf95_open("climoz_LMDZ.nc", nf90_nowrite, ncid)
    46 
    47        call nf95_inq_varid(ncid, "plev", varid)
    48        call nf95_gw_var(ncid, varid, plev)
    49        ! Convert from hPa to Pa because "paprs" and "pplay" are in Pa:
    50        plev = plev * 100.
    51        n_plev = size(plev)
    52     end if
    53 
    54     call bcast_mpi(n_plev)
    55     if (.not. is_mpi_root) allocate(plev(n_plev))
    56     call bcast_mpi(plev)
    57    
    58     ! Compute edges of pressure intervals:
    59     allocate(press_in_edg(n_plev + 1))
    60     if (is_mpi_root) then
    61        press_in_edg(1) = 0.
    62        ! We choose edges halfway in logarithm:
    63        forall (k = 2:n_plev) press_in_edg(k) = sqrt(plev(k - 1) * plev(k))
    64        press_in_edg(n_plev + 1) = huge(0.)
    65        ! (infinity, but any value guaranteed to be greater than the
    66        ! surface pressure would do)
    67     end if
    68     call bcast_mpi(press_in_edg)
    69     deallocate(plev) ! pointer
    70 
    71   end subroutine open_climoz
    72 
    73 end module open_climoz_m
     101END MODULE open_climoz_m

  • LMDZ5/branches/testing/libf/phylmd/phyaqua_mod.F90

    r2408 r2839  
    339339    PRINT *, 'iniaqua: before phyredem'
    340340
     341    pbl_tke(:,:,:)=1.e-8
    341342    falb1 = albedo
    342343    falb2 = albedo
    … …  
    354355    entr_therm = 0.
    355356    detr_therm = 0.
     357    ale_bl = 0.
     358    ale_bl_trig =0.
     359    alp_bl =0.
     360
    356361
    357362

  • LMDZ5/branches/testing/libf/phylmd/phys_cal_mod.F90

    r2542 r2839  
    99  INTEGER,SAVE :: day_cur       ! current day
    1010!$OMP THREADPRIVATE(day_cur)
    11   INTEGER,SAVE :: days_elapsed  ! number of whole days since start of the simulation
     11  INTEGER,SAVE :: days_elapsed  ! number of whole days since start of the current year
    1212!$OMP THREADPRIVATE(days_elapsed)
    1313  INTEGER,SAVE :: mth_len       ! number of days in the current month
    … …  
    3333!$OMP THREADPRIVATE(calend)
    3434
    35 
    3635CONTAINS
    3736 
    3837  SUBROUTINE phys_cal_init(annee_ref,day_ref)
    39   USE IOIPSL, ONLY:  ymds2ju
    40   USE ioipsl_getin_p_mod, ONLY: getin_p
    41   IMPLICIT NONE
     38
     39    USE IOIPSL, ONLY:  ymds2ju
     40    USE ioipsl_getin_p_mod, ONLY: getin_p
     41
     42    IMPLICIT NONE
    4243    INTEGER,INTENT(IN) :: annee_ref
    4344    INTEGER,INTENT(IN) :: day_ref
    … …  
    7980
    8081END MODULE phys_cal_mod
    81 

  • LMDZ5/branches/testing/libf/phylmd/phys_local_var_mod.F90

    r2787 r2839  
    163163      REAL, SAVE, ALLOCATABLE :: lcc3dstra(:,:)
    164164      !$OMP THREADPRIVATE(lcc3dstra)
     165      REAL, SAVE, ALLOCATABLE :: od443aer(:)
     166      !$OMP THREADPRIVATE(od443aer)
    165167      REAL, SAVE, ALLOCATABLE :: od550aer(:)
    166168      !$OMP THREADPRIVATE(od550aer)
    … …  
    207209      REAL, SAVE, ALLOCATABLE :: loaddust(:)
    208210      !$OMP THREADPRIVATE(loaddust)
     211      REAL, SAVE, ALLOCATABLE :: loadno3(:)
     212      !$OMP THREADPRIVATE(loadno3)
    209213      REAL, SAVE, ALLOCATABLE :: load_tmp1(:)
    210214      !$OMP THREADPRIVATE(load_tmp1)
    … …  
    213217      REAL, SAVE, ALLOCATABLE :: load_tmp3(:)
    214218      !$OMP THREADPRIVATE(load_tmp3)
    215       REAL, SAVE, ALLOCATABLE :: load_tmp4(:)
    216       !$OMP THREADPRIVATE(load_tmp4)
    217       REAL, SAVE, ALLOCATABLE :: load_tmp5(:)
    218       !$OMP THREADPRIVATE(load_tmp5)
    219       REAL, SAVE, ALLOCATABLE :: load_tmp6(:)
    220       !$OMP THREADPRIVATE(load_tmp6)
    221       REAL, SAVE, ALLOCATABLE :: load_tmp7(:)
    222       !$OMP THREADPRIVATE(load_tmp7)
    223219
    224220!IM ajout variables CFMIP2/CMIP5
    … …  
    352348!>jyg+nrlmd
    353349  !
    354       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: wbeff, zmax_th, zq2m, zt2m
    355 !$OMP THREADPRIVATE(wbeff, zmax_th, zq2m, zt2m)
     350      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: wbeff, convoccur, zmax_th, zq2m, zt2m
     351!$OMP THREADPRIVATE(wbeff, convoccur, zmax_th, zq2m, zt2m)
    356352      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zt2m_min_mon, zt2m_max_mon
    357353!$OMP THREADPRIVATE(zt2m_min_mon, zt2m_max_mon)
    … …  
    487483!$OMP THREADPRIVATE(budg_sed_part)
    488484#endif
     485      REAL, ALLOCATABLE, SAVE, DIMENSION(:) :: pr_tropopause
     486!$OMP THREADPRIVATE(pr_tropopause)
    489487
    490488CONTAINS
    … …  
    578576      allocate(lcc3dcon(klon, klev))
    579577      allocate(lcc3dstra(klon, klev))
     578      allocate(od443aer(klon))
    580579      allocate(od550aer(klon))
    581580      allocate(od865aer(klon))
    … …  
    600599      allocate(loadss(klon))
    601600      allocate(loaddust(klon))
     601      allocate(loadno3(klon))
    602602      allocate(load_tmp1(klon))
    603603      allocate(load_tmp2(klon))
    604604      allocate(load_tmp3(klon))
    605       allocate(load_tmp4(klon))
    606       allocate(load_tmp5(klon))
    607       allocate(load_tmp6(klon))
    608       allocate(load_tmp7(klon))
    609605
    610606!IM ajout variables CFMIP2/CMIP5
    … …  
    628624!!          Wake variables
    629625      ALLOCATE(ale_wake(klon), alp_wake(klon))
     626      ale_wake(:)=0.
    630627      ALLOCATE(wake_h(klon),wake_k(klon))
    631628      ALLOCATE(wake_omg(klon, klev))
    … …  
    677674      ALLOCATE(kh(klon), kh_x(klon), kh_w(klon))
    678675!
    679       ALLOCATE(wbeff(klon), zmax_th(klon))
     676      ALLOCATE(wbeff(klon), convoccur(klon), zmax_th(klon))
    680677      ALLOCATE(zq2m(klon), zt2m(klon), weak_inversion(klon))
    681678      ALLOCATE(zt2m_min_mon(klon), zt2m_max_mon(klon))
    … …  
    764761      ALLOCATE (vsed_aer(klon,klev))
    765762#endif
     763      ALLOCATE (pr_tropopause(klon))
    766764
    767765END SUBROUTINE phys_local_var_init
    … …  
    837835      deallocate(lcc3dcon)
    838836      deallocate(lcc3dstra)
     837      deallocate(od443aer)
    839838      deallocate(od550aer)
    840839      deallocate(od865aer)
    … …  
    859858      deallocate(loadss)
    860859      deallocate(loaddust)
     860      deallocate(loadno3)
    861861      deallocate(load_tmp1)
    862862      deallocate(load_tmp2)
    863863      deallocate(load_tmp3)
    864       deallocate(load_tmp4)
    865       deallocate(load_tmp5)
    866       deallocate(load_tmp6)
    867       deallocate(load_tmp7)
    868864      deallocate(du_gwd_hines,dv_gwd_hines,d_t_hin)
    869865      deallocate(d_q_ch4)
    … …  
    935931      DEALLOCATE(kh, kh_x, kh_w)
    936932!
    937       DEALLOCATE(wbeff, zmax_th)
     933      DEALLOCATE(wbeff, convoccur, zmax_th)
    938934      DEALLOCATE(zq2m, zt2m, weak_inversion)
    939935      DEALLOCATE(zt2m_min_mon, zt2m_max_mon)
    … …  
    10181014      DEALLOCATE (budg_sed_part)
    10191015#endif
     1016      DEALLOCATE (pr_tropopause)
    10201017
    10211018END SUBROUTINE phys_local_var_end

  • LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90

    r2787 r2839  
    669669  TYPE(ctrl_out), SAVE :: o_wbeff = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
    670670    'wbeff', 'Conv. updraft velocity at LFC (<100)', 'm/s', (/ ('', i=1, 10) /))
     671  TYPE(ctrl_out), SAVE :: o_convoccur = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     672    'convoccur', 'Convective occurence', '', (/ ('', i=1, 10) /))
    671673  TYPE(ctrl_out), SAVE :: o_prw = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11/), &
    672674    'prw', 'Precipitable water', 'kg/m2', (/ ('', i=1, 10) /))
    … …  
    10921094    'OD_10um_STRAT', 'Stratospheric Aerosol Optical depth at 10 um ', '1', (/ ('', i=1, 10) /))
    10931095!
     1096  TYPE(ctrl_out), SAVE :: o_od443aer = ctrl_out((/ 2, 6, 10, 10, 10, 10, 11, 11, 11, 11/), &
     1097    'od443aer', 'Total aerosol optical depth at 440nm', '-', (/ ('', i=1, 10) /))
    10941098  TYPE(ctrl_out), SAVE :: o_od550aer = ctrl_out((/ 2, 6, 10, 10, 10, 10, 11, 11, 11, 11/), &
    10951099    'od550aer', 'Total aerosol optical depth at 550nm', '-', (/ ('', i=1, 10) /))
    … …  
    11341138  TYPE(ctrl_out), SAVE :: o_loaddust = ctrl_out((/ 2, 6, 10, 10, 10, 10, 11, 11, 11, 11/), &
    11351139    'loaddust', 'Column Load of Dust ', 'kg/m2', (/ ('', i=1, 10) /))
     1140  TYPE(ctrl_out), SAVE :: o_loadno3 = ctrl_out((/ 2, 6, 10, 10, 10, 10, 11, 11, 11, 11/), &
     1141    'loadno3', 'Column Load of Nitrate ', 'kg/m2', (/ ('', i=1, 10) /))
    11361142  TYPE(ctrl_out), SAVE :: o_swtoaas_nat = ctrl_out((/ 4, 6, 10, 10, 10, 10, 11, 11, 11, 11/), &
    11371143    'swtoaas_nat', 'Natural aerosol radiative forcing all-sky at TOA', 'W/m2', (/ ('', i=1, 10) /))

  • LMDZ5/branches/testing/libf/phylmd/phys_output_mod.F90

    r2787 r2839  
    55  USE indice_sol_mod
    66  USE phys_output_var_mod
    7   USE aero_mod, only : naero_spc,name_aero
    87  USE phys_output_write_mod, ONLY : phys_output_write
    98  REAL, DIMENSION(nfiles),SAVE :: ecrit_files
    … …  
    4039    USE phys_cal_mod, only : hour, calend
    4140    USE mod_phys_lmdz_para
    42     USE aero_mod, only : naero_spc,name_aero
    4341    !Martin
    4442    USE surface_data, ONLY : ok_snow
    … …  
    4644    USE mod_grid_phy_lmdz, only: klon_glo,nbp_lon,nbp_lat
    4745    USE print_control_mod, ONLY: prt_level,lunout
    48     USE vertical_layers_mod, ONLY: ap,bp,preff,presnivs
     46    USE vertical_layers_mod, ONLY: ap,bp,preff,presnivs, aps, bps, pseudoalt
    4947    USE time_phylmdz_mod, ONLY: day_ini, itau_phy, start_time, annee_ref, day_ref
    5048#ifdef CPP_XIOS
    … …  
    9593    INTEGER                               :: idayref
    9694    REAL                                  :: zjulian_start, zjulian
    97     REAL, DIMENSION(klev)                 :: Ahyb, Bhyb, Alt
    9895    CHARACTER(LEN=4), DIMENSION(nlevSTD)  :: clevSTD
    9996    REAL, DIMENSION(nlevSTD)              :: rlevSTD
    … …  
    293290    zdtime_moy = dtime         ! Frequence ou l on moyenne
    294291
    295     ! Calcul des Ahyb, Bhyb et Alt
    296     DO k=1,klev
    297        Ahyb(k)=(ap(k)+ap(k+1))/2.
    298        Bhyb(k)=(bp(k)+bp(k+1))/2.
    299        Alt(k)=log(preff/presnivs(k))*8.
    300     ENDDO
    301     !          if(prt_level.ge.1) then
    302     WRITE(lunout,*)'Ap Hybrid = ',Ahyb(1:klev)
    303     WRITE(lunout,*)'Bp Hybrid = ',Bhyb(1:klev)
    304     WRITE(lunout,*)'Alt approx des couches pour une haut d echelle de 8km = ',Alt(1:klev)
    305     !          ENDIF
    306292
    307293  ecrit_files(7) = ecrit_files(1)
    … …  
    345331            levmax(iff) - levmin(iff) + 1, presnivs(levmin(iff):levmax(iff)))
    346332    CALL wxios_add_vaxis("Ahyb", &
    347             levmax(iff) - levmin(iff) + 1, Ahyb)
     333            levmax(iff) - levmin(iff) + 1, aps)
    348334    CALL wxios_add_vaxis("Bhyb", &
    349             levmax(iff) - levmin(iff) + 1, Bhyb)
     335            levmax(iff) - levmin(iff) + 1, bps)
    350336    CALL wxios_add_vaxis("Alt", &
    351             levmax(iff) - levmin(iff) + 1, Alt)
     337            levmax(iff) - levmin(iff) + 1, pseudoalt)
    352338   ELSE
    353339    ! NMC files
    … …  
    418404!!!! Composantes de la coordonnee sigma-hybride
    419405          CALL histvert(nid_files(iff), "Ahyb","Ahyb comp of Hyb Cord ", "Pa", &
    420                levmax(iff) - levmin(iff) + 1,Ahyb,nvertap(iff))
     406               levmax(iff) - levmin(iff) + 1,aps,nvertap(iff))
    421407
    422408          CALL histvert(nid_files(iff), "Bhyb","Bhyb comp of Hyb Cord", " ", &
    423                levmax(iff) - levmin(iff) + 1,Bhyb,nvertbp(iff))
     409               levmax(iff) - levmin(iff) + 1,bps,nvertbp(iff))
    424410
    425411          CALL histvert(nid_files(iff), "Alt","Height approx for scale heigh of 8km at levels", "Km", &                       
    426                levmax(iff) - levmin(iff) + 1,Alt,nvertAlt(iff))
     412               levmax(iff) - levmin(iff) + 1,pseudoalt,nvertAlt(iff))
    427413
    428414          ELSE

  • LMDZ5/branches/testing/libf/phylmd/phys_output_var_mod.F90

    r2787 r2839  
    2424  REAL, SAVE, ALLOCATABLE :: bils_latent(:) ! bilan de chaleur au sol
    2525  !$OMP THREADPRIVATE(bils_ec,bils_ech,bils_tke,bils_diss,bils_kinetic,bils_enthalp,bils_latent)
     26  ! output variables for energy conservation tests, computed in add_phys_tend
     27  REAL, SAVE, ALLOCATABLE :: d_qw_col(:)      ! watter vapour mass budget for each column (kg/m2/s)
     28  REAL, SAVE, ALLOCATABLE :: d_ql_col(:)      ! liquid watter mass budget for each column (kg/m2/s)
     29  REAL, SAVE, ALLOCATABLE :: d_qs_col(:)      ! solid watter mass budget for each column (kg/m2/s)
     30  REAL, SAVE, ALLOCATABLE :: d_qt_col(:)      ! total watter mass budget for each column (kg/m2/s)
     31  REAL, SAVE, ALLOCATABLE :: d_ek_col(:)      ! kinetic energy budget for each column (W/m2)
     32  REAL, SAVE, ALLOCATABLE :: d_h_dair_col(:)  ! enthalpy budget of dry air for each column (W/m2)
     33  REAL, SAVE, ALLOCATABLE :: d_h_qw_col(:)    ! enthalpy budget of watter vapour for each column (W/m2)
     34  REAL, SAVE, ALLOCATABLE :: d_h_ql_col(:)    ! enthalpy budget of liquid watter for each column (W/m2)
     35  REAL, SAVE, ALLOCATABLE :: d_h_qs_col(:)    ! enthalpy budget of solid watter  for each column (W/m2)
     36  REAL, SAVE, ALLOCATABLE :: d_h_col(:)       ! total enthalpy budget for each column (W/m2)
     37  !$OMP THREADPRIVATE(d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col)
     38  !$OMP THREADPRIVATE(d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col)
    2639
    2740! Marine
    … …  
    121134
    122135    allocate (bils_ec(klon),bils_ech(klon),bils_tke(klon),bils_diss(klon),bils_kinetic(klon),bils_enthalp(klon),bils_latent(klon))
     136    allocate (d_qw_col(klon), d_ql_col(klon), d_qs_col(klon), d_qt_col(klon), d_ek_col(klon), d_h_dair_col(klon) &
     137  &         , d_h_qw_col(klon), d_h_ql_col(klon), d_h_qs_col(klon), d_h_col(klon))
     138    d_qw_col=0. ; d_ql_col=0. ; d_qs_col=0. ; d_qt_col=0. ; d_ek_col=0. ; d_h_dair_col =0.
     139    d_h_qw_col=0. ; d_h_ql_col=0. ; d_h_qs_col=0. ; d_h_col=0.
    123140
    124141! Marine
    … …  
    155172    deallocate(snow_o,zfra_o,itau_con)
    156173    deallocate (bils_ec,bils_ech,bils_tke,bils_diss,bils_kinetic,bils_enthalp,bils_latent)
     174    deallocate (d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
     175  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col)
    157176
    158177! Marine

  • LMDZ5/branches/testing/libf/phylmd/phys_output_write_mod.F90

    r2787 r2839  
    6868         o_ue, o_ve, o_uq, o_vq, o_cape, o_pbase, &
    6969         o_ptop, o_fbase, o_plcl, o_plfc, &
    70          o_wbeff, o_cape_max, o_upwd, o_ep,o_epmax_diag, o_Ma, &
     70         o_wbeff, o_convoccur, o_cape_max, o_upwd, o_ep,o_epmax_diag, o_Ma, &
    7171         o_dnwd, o_dnwd0, o_ftime_con, o_mc, &
    7272         o_prw, o_prlw, o_prsw, o_s_pblh, o_s_pblt, o_s_lcl, &
    … …  
    9999         o_SWdownOR, o_LWdownOR, o_snowl, &
    100100         o_solldown, o_dtsvdfo, o_dtsvdft, &
    101          o_dtsvdfg, o_dtsvdfi, o_z0m, o_z0h, o_od550aer, &
     101         o_dtsvdfg, o_dtsvdfi, o_z0m, o_z0h, o_od443aer, o_od550aer, &
    102102         o_od865aer, o_absvisaer, o_od550lt1aer, &
    103103         o_sconcso4, o_sconcno3, o_sconcoa, o_sconcbc, &
    … …  
    105105         o_concoa, o_concbc, o_concss, o_concdust, &
    106106         o_loadso4, o_loadoa, o_loadbc, o_loadss, &
    107          o_loaddust, o_tausumaero, o_tausumaero_lw, &
     107         o_loaddust, o_loadno3, o_tausumaero, o_tausumaero_lw, &
    108108         o_topswad, o_topswad0, o_solswad, o_solswad0, &
    109109         o_toplwad, o_toplwad0, o_sollwad, o_sollwad0, &
    … …  
    196196#endif
    197197
    198     USE phys_state_var_mod, ONLY: pctsrf, paire_ter, rain_fall, snow_fall, &
     198    USE phys_state_var_mod, ONLY: pctsrf, rain_fall, snow_fall, &
    199199         qsol, z0m, z0h, fevap, agesno, &
    200200         nday_rain, rain_con, snow_con, &
    … …  
    236236         cldh, cldt, JrNt, cldljn, cldmjn, cldhjn, &
    237237         cldtjn, cldq, flwp, fiwp, ue, ve, uq, vq, &
    238          plcl, plfc, wbeff, upwd, dnwd, dnwd0, prw, prlw, prsw, &
     238         plcl, plfc, wbeff, convoccur, upwd, dnwd, dnwd0, prw, prlw, prsw, &
    239239         s_pblh, s_pblt, s_lcl, s_therm, uwriteSTD, &
    240240         vwriteSTD, wwriteSTD, phiwriteSTD, qwriteSTD, &
    … …  
    252252         weak_inversion, dthmin, cldtau, cldemi, &
    253253         pmflxr, pmflxs, prfl, psfl, re, fl, rh2m, &
    254          qsat2m, tpote, tpot, d_ts, od550aer, &
     254         qsat2m, tpote, tpot, d_ts, od443aer, od550aer, &
    255255         od865aer, absvisaer, od550lt1aer, sconcso4, sconcno3, &
    256256         sconcoa, sconcbc, sconcss, sconcdust, concso4, concno3, &
    257257         concoa, concbc, concss, concdust, loadso4, &
    258          loadoa, loadbc, loadss, loaddust, tausum_aero, &
     258         loadoa, loadbc, loadss, loaddust, loadno3, tausum_aero, &
    259259         topswad_aero, topswad0_aero, solswad_aero, &
    260260         solswad0_aero, topsw_aero, solsw_aero, &
    … …  
    325325    USE geometry_mod, ONLY: cell_area
    326326    USE surface_data, ONLY: type_ocean, version_ocean, ok_veget, ok_snow
    327 !    USE aero_mod, ONLY: naero_spc
    328327    USE aero_mod, ONLY: naero_tot, id_STRAT_phy
    329328    USE ioipsl, ONLY: histend, histsync
    … …  
    410409#ifdef CPP_XIOS
    411410#ifdef CPP_StratAer
    412     IF (.NOT.vars_defined) THEN
     411!$OMP MASTER
     412   IF (.NOT.vars_defined) THEN
    413413          !On ajoute les variables 3D traceurs par l interface fortran
    414414          CALL xios_get_handle("fields_strataer_trac_3D", group_handle)
    … …  
    474474          ENDDO
    475475    ENDIF
    476 #endif
    477 #endif
    478 
     476!$OMP END MASTER
     477#endif
     478#endif
    479479    ! ug la boucle qui suit ne sert qu'une fois, pour l'initialisation, sinon il n'y a toujours qu'un seul passage:
    480480    DO iinit=1, iinitend
    … …  
    513513       CALL histwrite_phy(o_contfracATM, zx_tmp_fi2d)
    514514       CALL histwrite_phy(o_contfracOR, pctsrf(:,is_ter))
    515        CALL histwrite_phy(o_aireTER, paire_ter)
    516515!
    517516#ifdef CPP_XIOS
    … …  
    739738       CALL histwrite_phy(o_bils_diss, bils_diss)
    740739       CALL histwrite_phy(o_bils_ec, bils_ec)
    741        IF (iflag_ener_conserv>=1) THEN
    742          CALL histwrite_phy(o_bils_ech, bils_ech)
    743        ENDIF
     740       CALL histwrite_phy(o_bils_ech, bils_ech)
    744741       CALL histwrite_phy(o_bils_tke, bils_tke)
    745742       CALL histwrite_phy(o_bils_kinetic, bils_kinetic)
    … …  
    890887             CALL histwrite_phy(o_plfc, plfc)
    891888             CALL histwrite_phy(o_wbeff, wbeff)
     889             CALL histwrite_phy(o_convoccur, convoccur)
    892890          ENDIF
    893891
    … …  
    11721170       IF (new_aod) THEN
    11731171          IF (flag_aerosol.GT.0) THEN
     1172             CALL histwrite_phy(o_od443aer, od443aer)
    11741173             CALL histwrite_phy(o_od550aer, od550aer)
    11751174             CALL histwrite_phy(o_od865aer, od865aer)
    … …  
    11931192             CALL histwrite_phy(o_loadss, loadss)
    11941193             CALL histwrite_phy(o_loaddust, loaddust)
     1194             CALL histwrite_phy(o_loadno3, loadno3)
    11951195             !--STRAT AER
    11961196          ENDIF

  • LMDZ5/branches/testing/libf/phylmd/phys_state_var_mod.F90

    r2787 r2839  
    275275      REAL,ALLOCATABLE,SAVE :: total_rain(:), nday_rain(:) 
    276276!$OMP THREADPRIVATE(total_rain,nday_rain)
    277       REAL,ALLOCATABLE,SAVE :: paire_ter(:)
    278 !$OMP THREADPRIVATE(paire_ter)
    279277! albsol1: albedo du sol total pour SW visible
    280278! albsol2: albedo du sol total pour SW proche IR
    … …  
    536534      ALLOCATE(pfrac_1nucl(klon,klev))
    537535      ALLOCATE(total_rain(klon), nday_rain(klon))
    538       ALLOCATE(paire_ter(klon))
    539536      ALLOCATE(albsol1(klon), albsol2(klon))
    540537!albedo SB >>>
    … …  
    675672      deallocate(pfrac_1nucl)
    676673      deallocate(total_rain, nday_rain)
    677       deallocate(paire_ter)
    678674      deallocate(albsol1, albsol2)
    679675!albedo SB >>>

  • LMDZ5/branches/testing/libf/phylmd/physiq_mod.F90

    r2787 r2839  
    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
    … …  
    154155!!!       d_s_the, d_dens_the, &            ! due to thermals
    155156       !                                 
    156        wbeff, zmax_th, &
     157       wbeff, convoccur, zmax_th, &
    157158       sens, flwp, fiwp,  &
    158159       ale_bl_stat,alp_bl_conv,alp_bl_det,  &
    … …  
    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
    … …  
    240242
    241243    USE cmp_seri_mod
     244    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, prt_enerbil, &
     245  &      fl_ebil, fl_cor_ebil
    242246
    243247    !IM stations CFMIP
    … …  
    245249    use FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
    246250    use ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
     251    USE VERTICAL_LAYERS_MOD, ONLY: aps,bps
     252
    247253
    248254    IMPLICIT none
    … …  
    404410    REAL pphis(klon)
    405411    REAL presnivs(klev)
    406     REAL znivsig(klev)
    407     real pir
     412!JLD    REAL znivsig(klev)
     413!JLD    real pir
    408414
    409415    REAL u(klon,klev)
    … …  
    465471    ! jmin_debut : indice minimum de j; nbptj : nombre de points en
    466472    ! direction j (latitude)
    467     INTEGER imin_debut, nbpti
    468     INTEGER jmin_debut, nbptj
     473!JLD    INTEGER imin_debut, nbpti
     474!JLD    INTEGER jmin_debut, nbptj
    469475    !IM: region='3d' <==> sorties en global
    470476    CHARACTER*3 region
    … …  
    530536    REAL Tconv(klon,klev)
    531537    REAL sij(klon,klev,klev)
     538!!    !
     539!!    ! variables pour tester la conservation de l'energie dans concvl
     540!!    REAL, DIMENSION(klon,klev)     :: d_t_con_sat
     541!!    REAL, DIMENSION(klon,klev)     :: d_q_con_sat
     542!!    REAL, DIMENSION(klon,klev)     :: dql_sat
    532543
    533544    real, save :: alp_bl_prescr=0.
    … …  
    615626    real env_tke_max0(klon)     ! TKE dans l'environnement au LCL
    616627
    617     !---D\'eclenchement stochastique
    618     integer :: tau_trig(klon)
     628!JLD    !---D\'eclenchement stochastique
     629!JLD    integer :: tau_trig(klon)
    619630
    620631    REAL,SAVE :: random_notrig_max=1.
    … …  
    650661    REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
    651662    ! RomP <<<
    652 
    653663    REAL          :: calday
    654664
    … …  
    749759    REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
    750760    !
     761#ifdef INCA
    751762    REAL zxsnow_dummy(klon)
     763#endif
    752764    REAL zsav_tsol(klon)
    753765    !
    … …  
    760772    LOGICAL zx_ajustq
    761773    !
    762     REAL za, zb
    763     REAL zx_t, zx_qs, zdelta, zcor, zlvdcp, zlsdcp
     774    REAL za
     775    REAL zx_t, zx_qs, zdelta, zcor
    764776    real zqsat(klon,klev)
    765777    !
    766     INTEGER i, k, iq, j, nsrf, ll, l
     778    INTEGER i, k, iq, nsrf, l
    767779    !
    768780    REAL t_coup
    … …  
    820832    REAL, dimension(klon)     :: dsig0, ddens0
    821833    INTEGER, dimension(klon)  :: wkoccur1
     834    ! tendance buffer pour appel de add_phys_tend
     835    REAL, DIMENSION(klon,klev)  :: d_q_ch4_dtime
    822836    !
    823837    ! Flag pour pouvoir ne pas ajouter les tendances.
    … …  
    872886
    873887    !
    874     integer itau_w   ! pas de temps ecriture = itap + itau_phy
     888!JLD    integer itau_w   ! pas de temps ecriture = itap + itau_phy
    875889    !
    876890    !
    … …  
    884898    REAL tabcntr0( length       )
    885899    !
    886     INTEGER ndex2d(nbp_lon*nbp_lat)
     900!JLD    INTEGER ndex2d(nbp_lon*nbp_lat)
    887901    !IM
    888902    !
    889903    !IM AMIP2 BEG
    890     REAL moyglo, mountor
     904!JLD    REAL moyglo, mountor
    891905    !IM 141004 BEG
    892906    REAL zustrdr(klon), zvstrdr(klon)
    … …  
    902916    !  REAL airedyn(nbp_lon+1,nbp_lat)
    903917    !IM 190504 END
    904     LOGICAL ok_msk
    905     REAL msk(klon)
     918!JLD    LOGICAL ok_msk
     919!JLD    REAL msk(klon)
    906920    !ym A voir plus tard
    907921    !ym      REAL zm_wo(jjmp1, klev)
    … …  
    910924    REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
    911925    REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
    912     REAL zx_tmp_2d(nbp_lon,nbp_lat)
    913     REAL zx_lon(nbp_lon,nbp_lat)
    914     REAL zx_lat(nbp_lon,nbp_lat)
     926!JLD    REAL zx_tmp_2d(nbp_lon,nbp_lat)
     927!JLD    REAL zx_lon(nbp_lon,nbp_lat)
     928!JLD    REAL zx_lat(nbp_lon,nbp_lat)
    915929    !
    916930    INTEGER nid_ctesGCM
    … …  
    928942    REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
    929943    !
    930     REAL zjulian
    931     SAVE zjulian
    932 !$OMP THREADPRIVATE(zjulian)
    933 
    934     INTEGER nhori, nvert
    935     REAL zsto
    936     REAL zstophy, zout
     944!JLD    REAL zjulian
     945!JLD    SAVE zjulian
     946!JLD!$OMP THREADPRIVATE(zjulian)
     947
     948!JLD    INTEGER nhori, nvert
     949!JLD    REAL zsto
     950!JLD    REAL zstophy, zout
    937951
    938952    character*20 modname
    … …  
    948962    parameter (prof2d_on = 1, prof3d_on = 2, &
    949963         prof2d_av = 3, prof3d_av = 4)
    950     !     Variables liees au bilan d'energie et d'enthalpi
    951964    REAL ztsol(klon)
    952     REAL      d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
    953     REAL      d_h_vcol_phy
    954     REAL      fs_bound, fq_bound
    955     SAVE      d_h_vcol_phy
    956     !$OMP THREADPRIVATE(d_h_vcol_phy)
    957     REAL      zero_v(klon)
    958     CHARACTER*40 ztit
    959     INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
    960     SAVE      ip_ebil
    961     DATA      ip_ebil/0/
    962     !$OMP THREADPRIVATE(ip_ebil)
    963965    REAL q2m(klon,nbsrf)  ! humidite a 2m
    964966
    965967    !IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
    966968    CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
    967     CHARACTER*40 tinst, tave, typeval
     969    CHARACTER*40 tinst, tave
    968970    REAL cldtaupi(klon,klev) ! Cloud optical thickness for
    969971    ! pre-industrial (pi) aerosols
    … …  
    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
    1027 
    1028     INTEGER ierr
     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.
     1026    ! vars_climoz(1:read_climoz-2) if read_climoz>2 (temporary)
     1027    REAL :: ro3i ! 0<=ro3i<=360 ; required time index in NetCDF file for
     1028                 ! the ozone fields, old method.
     1029
    10291030    include "YOMCST.h"
    10301031    include "YOETHF.h"
    … …  
    10401041    ! Declarations pour Simulateur COSP
    10411042    !============================================================
     1043#ifdef CPP_COSP
    10421044    real :: mr_ozone(klon,klev)
    1043 
     1045#endif
    10441046    !IM stations CFMIP
    10451047    INTEGER, SAVE :: nCFMIP
    … …  
    10901092    !--OB variables for mass fixer (hard coded for now)
    10911093    logical, parameter :: mass_fixer=.false.
    1092     real qql1(klon),qql2(klon),zdz,corrqql
     1094    real qql1(klon),qql2(klon),corrqql
    10931095
    10941096    ! Ehouarn: set value of jjmp1 since it is no longer a "fixed parameter"
    … …  
    11941196
    11951197    modname = 'physiq'
    1196     !IM
    1197     IF (ip_ebil_phy.ge.1) THEN
    1198        DO i=1,klon
    1199           zero_v(i)=0.
    1200        ENDDO
    1201     ENDIF
    12021198
    12031199    IF (debut) THEN
    … …  
    12111207       iflag_wake_tend = 0
    12121208       CALL getin_p('iflag_wake_tend',iflag_wake_tend)
     1209       ok_bad_ecmwf_thermo=.TRUE. ! By default thermodynamical constants are set
     1210                                  ! in rrtm/suphec.F90 (and rvtmp2 is set to 0).
     1211       CALL getin_p('ok_bad_ecmwf_thermo',ok_bad_ecmwf_thermo)
     1212       fl_ebil = 0 ! by default, conservation diagnostics are desactivated
     1213       CALL getin_p('fl_ebil',fl_ebil)
     1214       fl_cor_ebil = 0 ! by default, no correction to ensure energy conservation
     1215       CALL getin_p('fl_cor_ebil',fl_cor_ebil)
    12131216    ENDIF
    12141217
    … …  
    12741277       clwcon(:,:) = 0.0
    12751278
    1276        !IM     
    1277        IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
    12781279       !
    12791280       print*,'iflag_coupl,iflag_clos,iflag_wake', &
    … …  
    16561657               start_time, &
    16571658               itau_phy, &
     1659               date0, &
    16581660               io_lon, &
    16591661               io_lat)
    … …  
    16711673
    16721674       !$omp single
    1673        IF (read_climoz >= 1) THEN
    1674           CALL open_climoz(ncid_climoz, press_climoz)
    1675        ENDIF
     1675       IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz,   &
     1676           press_edg_climoz, time_climoz, ok_daily_climoz, adjust_tropopause)
    16761677       !$omp end single
    16771678       !
    … …  
    17371738    !
    17381739    itap   = itap + 1
    1739     IF (is_mpi_root .AND. is_omp_root) THEN
     1740    IF (is_master .OR. prt_level > 9) THEN
    17401741      IF (prt_level > 5 .or. MOD(itap,5) == 0) THEN
    1741          WRITE(LUNOUT,*)'Entering physics current time = ', current_time
    1742          WRITE(LUNOUT,*)'Date = ',year_cur,'/',mth_cur,'/',day_cur,':',hour/3600
     1742         WRITE(LUNOUT,*)'Entering physics elapsed seconds since start ', current_time
     1743         WRITE(LUNOUT,100)year_cur,mth_cur,day_cur,hour/3600.
     1744 100     FORMAT('Date = ',i4.4,' / ',i2.2, ' / ',i2.2,' : ',f20.17)
    17431745      ENDIF
    17441746    ENDIF
    … …  
    19781980          wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
    19791981       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)
     1982          !--- ro3i = elapsed days number since current year 1st january, 0h
     1983          ro3i=days_elapsed+jh_cur-jh_1jan
     1984          !--- scaling for old style files (360 records)
     1985          IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
     1986          IF(adjust_tropopause) THEN
     1987             CALL dyn_tropopause(t_seri, ztsol, paprs, pplay, presnivs, rot, &
     1988                  pr_tropopause)
     1989             CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),   &
     1990                      ro3i,  press_edg_climoz,  paprs,   wo, time_climoz,    &
     1991                           latitude_deg,  press_cen_climoz,  pr_tropopause)
    19891992          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
     1993             CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),   &
     1994                      ro3i,  press_edg_climoz,  paprs,  wo,  time_climoz)
     1995          END IF
    19951996          ! Convert from mole fraction of ozone to column density of ozone in a
    19961997          ! cell, in kDU:
    19971998          FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
    19981999               * zmasse / dobson_u / 1e3
    1999           ! (By regridding ozone values for LMDZ only once per day, we
     2000          ! (By regridding ozone values for LMDZ only once a day, we
    20002001          ! have already neglected the variation of pressure in one
    20012002          ! day. So do not recompute "wo" at each time step even if
    … …  
    20112012     CALL add_phys_tend &
    20122013            (du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
    2013                'eva',abortphy,flag_inhib_tend)
     2014               'eva',abortphy,flag_inhib_tend,itap,0)
     2015    call prt_enerbil('eva',itap)
    20142016
    20152017    !=========================================================================
    … …  
    22332235          CALL add_pbl_tend &
    22342236               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
    2235                'vdf',abortphy,flag_inhib_tend)
     2237               'vdf',abortphy,flag_inhib_tend,itap)
    22362238       ELSE
    22372239          CALL add_phys_tend &
    22382240               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
    2239                'vdf',abortphy,flag_inhib_tend)
    2240        ENDIF
     2241               'vdf',abortphy,flag_inhib_tend,itap,0)
     2242       ENDIF
     2243       call prt_enerbil('vdf',itap)
    22412244       !--------------------------------------------------------------------
    22422245
    … …  
    24852488               d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
    24862489               rain_con, snow_con, ibas_con, itop_con, sigd, &
    2487                ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0, &
     2490               ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
    24882491               Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
    24892492               pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
    … …  
    26182621    itapcv = itapcv+1
    26192622
     2623!!!jyg  Appel diagnostique a add_phys_tend pour tester la conservation de
     2624!!!     l'energie dans les courants satures.
     2625!!    d_t_con_sat(:,:) = d_t_con(:,:) - ftd(:,:)*dtime
     2626!!    d_q_con_sat(:,:) = d_q_con(:,:) - fqd(:,:)*dtime
     2627!!    dql_sat(:,:) = (wdtrainA(:,:)+wdtrainM(:,:))*dtime/zmasse(:,:)
     2628!!    CALL add_phys_tend(d_u_con, d_v_con, d_t_con_sat, d_q_con_sat, dql_sat,   &
     2629!!                     dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,& 
     2630!!                     itap, 1)
     2631!!    call prt_enerbil('convection_sat',itap)
     2632!!
     2633!!
    26202634    CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
    2621          'convection',abortphy,flag_inhib_tend)
     2635         'convection',abortphy,flag_inhib_tend,itap,0)
     2636    call prt_enerbil('convection',itap)
    26222637
    26232638    !-------------------------------------------------------------------------
    … …  
    27482763       ! ajout des tendances des poches froides
    27492764       CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
    2750             abortphy,flag_inhib_tend)
     2765            abortphy,flag_inhib_tend,itap,0)
     2766       call prt_enerbil('wake',itap)
    27512767       !------------------------------------------------------------------------
    27522768
    … …  
    27572773            (d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk, wake_k, &
    27582774             'wake', abortphy)
    2759 
     2775          call prt_enerbil('wake',itap)
    27602776       ENDIF   ! (iflag_wake_tend .GT. 0.)
    27612777
    … …  
    28782894             CALL add_wake_tend &
    28792895                 (d_deltat_the, d_deltaq_the, dsig0, ddens0, wkoccur1, 'the', abortphy)
     2896             call prt_enerbil('the',itap)
    28802897          !
    28812898          ENDIF  ! (mod(iflag_pbl_split/2,2) .EQ. 1)
    28822899          !
    28832900          CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,  &
    2884                              dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend)
     2901                             dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
     2902          call prt_enerbil('thermals',itap)
    28852903          !
    28862904!
    … …  
    29442962          ! ajout des tendances de l'ajustement sec ou des thermiques
    29452963          CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
    2946                'ajsb',abortphy,flag_inhib_tend)
     2964               'ajsb',abortphy,flag_inhib_tend,itap,0)
     2965          call prt_enerbil('ajsb',itap)
    29472966          d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
    29482967          d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
    … …  
    29873006    WHERE (snow_lsc < 0) snow_lsc = 0.
    29883007
     3008!+JLD
     3009!    write(*,9000) 'phys lsc',"enerbil: bil_q, bil_e,",rain_lsc+snow_lsc &
     3010!        & ,((rcw-rcpd)*rain_lsc + (rcs-rcpd)*snow_lsc)*t_seri(1,1)-rlvtt*rain_lsc+rlstt*snow_lsc &
     3011!        & ,rain_lsc,snow_lsc
     3012!    write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
     3013!-JLD
    29893014    CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
    2990          'lsc',abortphy,flag_inhib_tend)
     3015         'lsc',abortphy,flag_inhib_tend,itap,0)
     3016    call prt_enerbil('lsc',itap)
    29913017    rain_num(:)=0.
    29923018    DO k = 1, klev
    … …  
    37673793    ENDDO
    37683794
    3769     CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend)
    3770     CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend)
     3795    CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
     3796    call prt_enerbil('SW',itap)
     3797    CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
     3798    call prt_enerbil('LW',itap)
    37713799
    37723800    !
    … …  
    38383866       ! ajout des tendances de la trainee de l'orographie
    38393867       CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
    3840             abortphy,flag_inhib_tend)
     3868            abortphy,flag_inhib_tend,itap,0)
     3869       call prt_enerbil('oro',itap)
    38413870       !----------------------------------------------------------------------
    38423871       !
    … …  
    38843913       ! ajout des tendances de la portance de l'orographie
    38853914       CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
    3886             'lif', abortphy,flag_inhib_tend)
     3915            'lif', abortphy,flag_inhib_tend,itap,0)
     3916       call prt_enerbil('lif',itap)
    38873917    ENDIF ! fin de test sur ok_orolf
    38883918
    … …  
    39073937       d_t_hin(:, :)=0.
    39083938       CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
    3909             dqi0, paprs, 'hin', abortphy,flag_inhib_tend)
     3939            dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
     3940       call prt_enerbil('hin',itap)
    39103941    ENDIF
    39113942
    … …  
    39243955
    39253956       CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
    3926             paprs, 'front_gwd_rando', abortphy,flag_inhib_tend)
     3957            paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
     3958       call prt_enerbil('front_gwd_rando',itap)
    39273959    ENDIF
    39283960
    … …  
    39323964            du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
    39333965       CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
    3934             paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend)
     3966            paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
     3967       call prt_enerbil('flott_gwd_rando',itap)
    39353968       zustr_gwd_rando=0.
    39363969       zvstr_gwd_rando=0.
    … …  
    39814014       CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
    39824015       ! ajout de la tendance d'humidite due au methane
    3983        CALL add_phys_tend(du0, dv0, dt0, d_q_ch4*dtime, dql0, dqi0, paprs, &
    3984             'q_ch4', abortphy,flag_inhib_tend)
     4016       d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*dtime
     4017       CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
     4018            'q_ch4', abortphy,flag_inhib_tend,itap,0)
     4019       d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/dtime
    39854020    ENDIF
    39864021    !
    … …  
    40044039          CALL phys_cosp(itap,dtime,freq_cosp, &
    40054040               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
    4006                ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, &
     4041               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
    40074042               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
    40084043               JrNt,ref_liq,ref_ice, &
    … …  
    42114246            pphi, &
    42124247            pphis, &
    4213             zx_rh)
     4248            zx_rh, &
     4249            aps, bps)
    42144250
    42154251       CALL VTe(VTinca)
    … …  
    44654501             CALL nf95_close(ncid_climoz)
    44664502          ENDIF
    4467           DEALLOCATE(press_climoz) ! pointer
     4503          DEALLOCATE(press_edg_climoz) ! pointer
     4504          DEALLOCATE(press_cen_climoz) ! pointer
    44684505       ENDIF
    44694506       !$OMP END MASTER

  • LMDZ5/branches/testing/libf/phylmd/readaerosol.F90

    r2408 r2839  
    222222    REAL, ALLOCATABLE, DIMENSION(:)       :: varktmp
    223223
    224     REAL, DIMENSION(nbp_lon,nbp_lat,12)         :: psurf_glo2D   ! Surface pression for 12 months on dynamics global grid
     224    REAL, DIMENSION(nbp_lon,nbp_lat,12)   :: psurf_glo2D   ! Surface pression for 12 months on dynamics global grid
    225225    REAL, DIMENSION(klon_glo,12)          :: psurf_glo1D   ! -"- on physical global grid
    226     REAL, DIMENSION(nbp_lon,nbp_lat,12)         :: load_glo2D    ! Load for 12 months on dynamics global grid
     226    REAL, DIMENSION(nbp_lon,nbp_lat,12)   :: load_glo2D    ! Load for 12 months on dynamics global grid
    227227    REAL, DIMENSION(klon_glo,12)          :: load_glo1D    ! -"- on physical global grid
    228     REAL, DIMENSION(nbp_lon,nbp_lat)            :: vartmp
    229     REAL, DIMENSION(nbp_lon)                  :: lon_src              ! longitudes in file
    230     REAL, DIMENSION(nbp_lat)                :: lat_src, lat_src_inv ! latitudes in file
     228    REAL, DIMENSION(nbp_lon,nbp_lat)      :: vartmp
     229    REAL, DIMENSION(nbp_lon)              :: lon_src              ! longitudes in file
     230    REAL, DIMENSION(nbp_lat)              :: lat_src, lat_src_inv ! latitudes in file
    231231    LOGICAL                               :: new_file             ! true if new file format detected
    232232    LOGICAL                               :: invert_lat           ! true if the field has to be inverted for latitudes
    … …  
    245245 
    246246       WRITE(lunout,*) 'reading variable ',TRIM(varname),' in file ', TRIM(fname)
    247        CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid), "pb open "//trim(varname) )
     247       CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid), "pb open "//trim(fname) )
    248248
    249249! Test for equal longitudes and latitudes in file and model
    … …  
    338338!****************************************************************************************
    339339          ! Get variable id
    340           CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid),"pb inq var "//TRIM(varname) )
    341          
    342           ! Get the variable
    343           CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)),"pb get var "//TRIM(varname) )
     340          !CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid),"pb inq var "//TRIM(varname) )
     341          print *,'readaerosol ', TRIM(varname)
     342          IF ( nf90_inq_varid(ncid, TRIM(varname), varid) /= NF90_NOERR ) THEN
     343            ! Variable is not there
     344            WRITE(lunout,*) 'Attention '//TRIM(varname)//' is not in aerosol input file'
     345            varyear(:,:,:,:)=0.0
     346          ELSE
     347            ! Get the variable
     348            CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)),"pb get var "//TRIM(varname) )
     349          ENDIF
    344350         
    345351! ++) Read surface pression, 12 month in one variable
    … …  
    353359!****************************************************************************************
    354360          ! Get variable id
    355           CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) ,"pb inq var load_"//TRIM(varname))
    356           ! Get the variable
    357           CALL check_err( nf90_get_var(ncid, varid, load_glo2D),"pb get var load_"//TRIM(varname) )
     361          !CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) ,"pb inq var load_"//TRIM(varname))
     362          IF ( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) /= NF90_NOERR ) THEN
     363            WRITE(lunout,*) 'Attention load_'//TRIM(varname)//' is not in aerosol input file'
     364            load_glo2D(:,:,:)=0.0
     365          ELSE
     366            ! Get the variable
     367            CALL check_err( nf90_get_var(ncid, varid, load_glo2D),"pb get var load_"//TRIM(varname) )
     368          ENDIF
    358369         
    359370! ++) Read ap

  • LMDZ5/branches/testing/libf/phylmd/readaerosol_optic.F90

    r2669 r2839  
    1616  USE phys_local_var_mod, only: sconcso4,sconcno3,sconcoa,sconcbc,sconcss,sconcdust, &
    1717      concso4,concno3,concoa,concbc,concss,concdust,loadso4,loadoa,loadbc,loadss,loaddust, &
    18       load_tmp1,load_tmp2,load_tmp3,load_tmp4,load_tmp5,load_tmp6,load_tmp7
     18      load_tmp1,load_tmp2,load_tmp3
    1919  IMPLICIT NONE
    2020
    … …  
    9393
    9494     ! Get bc aerosol distribution
    95      CALL readaerosol_interp(id_ASBCM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcsol, bcsol_pi, load_tmp1 )
    96      CALL readaerosol_interp(id_AIBCM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcins, bcins_pi, load_tmp2 )
     95     CALL readaerosol_interp(id_ASBCM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcsol, bcsol_pi, load_tmp1)
     96     CALL readaerosol_interp(id_AIBCM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcins, bcins_pi, load_tmp2)
    9797     loadbc(:)=load_tmp1(:)+load_tmp2(:)
    9898  ELSE
    … …  
    107107       flag_aerosol .EQ. 6 ) THEN
    108108
    109      CALL readaerosol_interp(id_ASPOMM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomsol, pomsol_pi, load_tmp3)
    110      CALL readaerosol_interp(id_AIPOMM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomins, pomins_pi, load_tmp4)
    111      loadoa(:)=load_tmp3(:)+load_tmp4(:)
     109     CALL readaerosol_interp(id_ASPOMM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomsol, pomsol_pi, load_tmp1)
     110     CALL readaerosol_interp(id_AIPOMM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomins, pomins_pi, load_tmp2)
     111     loadoa(:)=load_tmp1(:)+load_tmp2(:)
    112112  ELSE
    113113     pomsol(:,:) = 0. ; pomsol_pi(:,:) = 0.
    … …  
    121121      flag_aerosol .EQ. 6 ) THEN
    122122
    123       CALL readaerosol_interp(id_SSSSM_phy ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sssupco, sssupco_pi, load_tmp5)
    124       CALL readaerosol_interp(id_CSSSM_phy ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sscoarse,sscoarse_pi, load_tmp6)
    125       CALL readaerosol_interp(id_ASSSM_phy ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, ssacu, ssacu_pi, load_tmp7)
    126      loadss(:)=load_tmp5(:)+load_tmp6(:)+load_tmp7(:)
    127   ELSE
    128      sscoarse(:,:) = 0. ; sscoarse_pi(:,:) = 0.
    129      ssacu(:,:)    = 0. ; ssacu_pi(:,:) = 0.
    130      sssupco(:,:)  = 0. ; sssupco_pi = 0.
    131      loadss=0.
     123      CALL readaerosol_interp(id_SSSSM_phy ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sssupco, sssupco_pi,  load_tmp1)
     124      CALL readaerosol_interp(id_CSSSM_phy ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sscoarse,sscoarse_pi, load_tmp2)
     125      CALL readaerosol_interp(id_ASSSM_phy ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, ssacu,   ssacu_pi,    load_tmp3)
     126      loadss(:)=load_tmp1(:)+load_tmp2(:)+load_tmp3(:)
     127  ELSE
     128      sscoarse(:,:) = 0. ; sscoarse_pi(:,:) = 0.
     129      ssacu(:,:)    = 0. ; ssacu_pi(:,:) = 0.
     130      sssupco(:,:)  = 0. ; sssupco_pi = 0.
     131      loadss=0.
    132132  ENDIF
    133133
    … …  
    137137
    138138      CALL readaerosol_interp(id_CIDUSTM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, cidust, cidust_pi, loaddust)
    139 
    140139  ELSE
    141140      cidust(:,:) = 0. ; cidust_pi(:,:) = 0.
    … …  
    174173  m_allaer_pi(:,:,id_CSNO3M_phy) = 0.0
    175174  m_allaer_pi(:,:,id_CINO3M_phy) = 0.0
    176 
    177175!
    178176! Calculate the total mass of all soluble aersosols
    … …  
    214212     
    215213  END IF
    216 
    217214
    218215! Diagnostics calculation for CMIP5 protocol
    … …  
    230227  concdust(:,:)=m_allaer(:,:,id_CIDUSTM_phy)*1.e-9
    231228
    232 
    233229END SUBROUTINE readaerosol_optic

  • LMDZ5/branches/testing/libf/phylmd/reevap.F90

    r2720 r2839  
    22   &         d_t_eva,d_q_eva,d_ql_eva,d_qs_eva)
    33
    4 
     4    ! flag to include modifications to ensure energy conservation (if flag >0)
     5    USE add_phys_tend_mod, only : fl_cor_ebil
     6   
    57    IMPLICIT none
    68    !>======================================================================
    … …  
    2224    ! Re-evaporer l'eau liquide nuageuse
    2325    !
    24 
     26!print *,'rrevap ; fl_cor_ebil:',fl_cor_ebil,' iflag_ice_thermo:',iflag_ice_thermo,' RVTMP2',RVTMP2
    2527    DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
    2628       DO i = 1, klon
     29        if (fl_cor_ebil .GT. 0) then
     30          zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*(q_seri(i,k)+ql_seri(i,k)+qs_seri(i,k)))
     31          zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*(q_seri(i,k)+ql_seri(i,k)+qs_seri(i,k)))
     32        else
    2733          zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
    2834          !jyg<
    … …  
    3036          !                  A verifier !!!
    3137          zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
     38        end if
    3239          IF (iflag_ice_thermo .EQ. 0) THEN
    3340             zlsdcp=zlvdcp
    … …  
    3946
    4047             zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
     48  zdelta = 0.
    4149             zb = MAX(0.0,ql_seri(i,k))
    4250             za = - MAX(0.0,ql_seri(i,k)) &

  • LMDZ5/branches/testing/libf/phylmd/regr_lat_time_coefoz_m.F90

    r2471 r2839  
    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/branches/testing/libf/phylmd/regr_pr_comb_coefoz_m.F90

    r2408 r2839  
    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/branches/testing/libf/phylmd/regr_pr_int_m.F90

    r2408 r2839  
    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/branches/testing/libf/phylmd/regr_pr_o3_m.F90

    r2471 r2839  
    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)

  • LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90

    r2787 r2839  
    1212  USE DIMPHY
    1313  USE aero_mod
    14   USE phys_local_var_mod, ONLY: od550aer,od865aer,ec550aer,od550lt1aer
     14  USE phys_local_var_mod, ONLY: od443aer,od550aer,od865aer,ec550aer,od550lt1aer
    1515  USE YOMCST, ONLY: RD,RG
    1616
    … …  
    327327        soluble=.TRUE.
    328328        spsol=4
    329         fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD
     329        !fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD
     330        fac=0.0      !--6 March 2017 - OB as Didier H said CSSO4 should not be used
    330331    ELSEIF (aerosol_name(m).EQ.id_SSSSM_phy) THEN
    331332        soluble=.TRUE.
    … …  
    366367    DO la=1,las
    367368
    368     !--only 550, 670 and 865 nm are used
    369     IF (la.NE.la550.AND.la.NE.la670.AND.la.NE.la865) CYCLE
     369    !--only 443, 550, 670 and 865 nm are used
     370    !--to save time 670 and AI are not computed for CMIP6
     371    !IF (la.NE.la443.AND.la.NE.la550.AND.la.NE.la670.AND.la.NE.la865) CYCLE
     372    IF (la.NE.la443.AND.la.NE.la550.AND.la.NE.la865) CYCLE
    370373
    371374      IF (soluble) THEN            ! For soluble aerosol
    … …  
    433436
    434437!--AOD calculations for diagnostics
     438  od443aer(:)=SUM(tausum(:,la443,:),dim=2)
    435439  od550aer(:)=SUM(tausum(:,la550,:),dim=2)
    436   od670aer(:)=SUM(tausum(:,la670,:),dim=2)
     440  !od670aer(:)=SUM(tausum(:,la670,:),dim=2)
    437441  od865aer(:)=SUM(tausum(:,la865,:),dim=2)
    438442
    … …  
    441445
    442446!--aerosol index
    443   ai(:)=-LOG(MAX(od670aer(:),1.e-8)/MAX(od865aer(:),1.e-8))/LOG(670./865.)
     447  ai(:)=0.0
     448  !ai(:)=-LOG(MAX(od670aer(:),1.e-8)/MAX(od865aer(:),1.e-8))/LOG(670./865.)
    444449
    445450  od550lt1aer(:)=tausum(:,la550,id_ASSO4M_phy)+tausum(:,la550,id_ASBCM_phy) +tausum(:,la550,id_AIBCM_phy)+ &

  • LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90

    r2787 r2839  
    566566        soluble=.TRUE.
    567567        spsol=3
    568         fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD
     568        !fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD
     569        fac=0.0      !--6 March 2017 - OB as Didier H said CSSO4 should not be used
    569570     ELSEIF  (aerosol_name(m).EQ.id_ASSO4M_phy) THEN
    570571        soluble=.TRUE.

  • LMDZ5/branches/testing/libf/phylmd/rrtm/read_rsun_rrtm.F90

    r2787 r2839  
    1010  USE netcdf, ONLY: nf90_get_var, nf90_noerr, nf90_nowrite
    1111
    12   USE phys_cal_mod, ONLY : day_cur, year_len
     12  USE phys_cal_mod, ONLY : days_elapsed, year_len
    1313
    1414  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
    … …  
    2323
    2424  ! Input arguments
    25   LOGICAL, INTENT(IN)          :: debut
     25  LOGICAL, INTENT(IN) :: debut
    2626
    2727! Local variables
    28   INTEGER               :: ncid, dimid, varid, ncerr, nbday
     28  INTEGER :: ncid, dimid, varid, ncerr, nbday
    2929  REAL, POINTER :: wlen(:), time(:)
    30   REAL, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: SSI_FRAC
     30  REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: SSI_FRAC
    3131!$OMP THREADPRIVATE(SSI_FRAC)
    3232  REAL, ALLOCATABLE, SAVE, DIMENSION(:) :: TSI(:)
    … …  
    9090    ENDIF
    9191
    92 !--only read at beginning of month
    93     IF (debut.OR.day_cur.NE.day_pre) THEN
     92!--only read at beginning of day
     93!--day in year is provided as days_elapsed since the beginning of the year +1
     94    IF (debut.OR.days_elapsed+1.NE.day_pre) THEN
    9495
    95 !--keep memory of previous month
    96       day_pre=day_cur
     96!--keep memory of previous day
     97      day_pre=days_elapsed+1
    9798
    9899!--copy
    99       RSUN(1:NSW)=SSI_FRAC(:,day_cur)
    100       solaire=TSI(day_cur)
     100      RSUN(1:NSW)=SSI_FRAC(:,days_elapsed+1)
     101      solaire=TSI(days_elapsed+1)
    101102
    102       print *,'READ_RSUN_RRTM day=', day_cur,' solaire=', solaire, ' RSUN=', RSUN(1:NSW)
     103      print *,'READ_RSUN_RRTM day=', days_elapsed+1,' solaire=', solaire, ' RSUN=', RSUN(1:NSW)
    103104
    104105    ENDIF !--fin allocation

  • LMDZ5/branches/testing/libf/phylmd/rrtm/readaerosol_optic_rrtm.F90

    r2787 r2839  
    1717  USE phys_local_var_mod, only: sconcso4,sconcno3,sconcoa,sconcbc,sconcss,sconcdust, &
    1818       concso4,concno3,concoa,concbc,concss,concdust,loadso4,loadoa,loadbc,loadss,loaddust, &
    19        load_tmp1,load_tmp2,load_tmp3,load_tmp4,load_tmp5,load_tmp6,load_tmp7
     19       loadno3, load_tmp1,load_tmp2,load_tmp3
    2020
    2121  USE infotrac_phy
    … …  
    195195     IF ( flag_aerosol .EQ. 3 .OR. flag_aerosol .EQ. 6 ) THEN
    196196
    197         CALL readaerosol_interp(id_ASPOMM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomsol, pomsol_pi, load_tmp3)
    198         CALL readaerosol_interp(id_AIPOMM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomins, pomins_pi, load_tmp4)
    199         loadoa(:)=load_tmp3(:)+load_tmp4(:)
     197        CALL readaerosol_interp(id_ASPOMM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomsol, pomsol_pi, load_tmp1)
     198        CALL readaerosol_interp(id_AIPOMM_phy, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomins, pomins_pi, load_tmp2)
     199        loadoa(:)=load_tmp1(:)+load_tmp2(:)
    200200     ELSE
    201201        pomsol(:,:) = 0. ; pomsol_pi(:,:) = 0.
    … …  
    208208
    209209        CALL readaerosol_interp(id_SSSSM_phy ,itap, pdtphys,rjourvrai, &
    210         debut, pplay, paprs, t_seri, sssupco, sssupco_pi, load_tmp5)
     210        debut, pplay, paprs, t_seri, sssupco, sssupco_pi, load_tmp1)
    211211        CALL readaerosol_interp(id_CSSSM_phy ,itap, pdtphys,rjourvrai, &
    212         debut, pplay, paprs, t_seri, sscoarse,sscoarse_pi, load_tmp6)
    213         CALL readaerosol_interp(id_ASSSM_phy ,itap, pdtphys, rjourvrai, &
    214         debut, pplay, paprs, t_seri, ssacu, ssacu_pi, load_tmp7)
    215         loadss(:)=load_tmp5(:)+load_tmp6(:)+load_tmp7(:)
     212        debut, pplay, paprs, t_seri, sscoarse,sscoarse_pi, load_tmp2)
     213        CALL readaerosol_interp(id_ASSSM_phy ,itap, pdtphys,rjourvrai, &
     214        debut, pplay, paprs, t_seri, ssacu, ssacu_pi, load_tmp3)
     215        loadss(:)=load_tmp1(:)+load_tmp2(:)+load_tmp3(:)
    216216     ELSE
    217217        sscoarse(:,:) = 0. ; sscoarse_pi(:,:) = 0.
    … …  
    231231     ENDIF
    232232     !
     233     ! Read and interpolate cidustm
     234     IF (flag_aerosol .EQ. 6) THEN
     235
     236        CALL readaerosol_interp(id_ASNO3M_phy, itap, pdtphys, rjourvrai, &
     237        debut, pplay, paprs, t_seri, nitracc, nitracc_pi, load_tmp1)
     238        CALL readaerosol_interp(id_CSNO3M_phy, itap, pdtphys, rjourvrai, &
     239        debut, pplay, paprs, t_seri, nitrcoarse, nitrcoarse_pi, load_tmp2)
     240        CALL readaerosol_interp(id_CINO3M_phy, itap, pdtphys, rjourvrai, &
     241        debut, pplay, paprs, t_seri, nitrinscoarse, nitrinscoarse_pi, load_tmp3)
     242        loadss(:)=load_tmp1(:)+load_tmp2(:)+load_tmp3(:)
     243
     244     ELSE
     245        nitracc(:,:)         =   0.0 ; nitracc_pi(:,:)      =   0.0
     246        nitrcoarse(:,:)      =   0.0 ; nitrcoarse_pi(:,:)   =   0.0
     247        nitrinscoarse(:,:)   =   0.0 ; nitrinscoarse_pi(:,:)=   0.0
     248        loadno3(:)=0.0
     249     ENDIF
     250     !
     251     ! CSSO4M is set to 0 as not reliable
    233252     sulfcoarse(:,:)      =   0.0 ! CSSO4M (=SO4) + CSMSAM (=MSA)
    234253     sulfcoarse_pi(:,:)   =   0.0 ! CSSO4M (=SO4) + CSMSAM (=MSA) pre-ind
    235      !
    236      !--placeholder for offline nitrate   
    237      !
    238      nitracc(:,:)         =   0.0
    239      nitracc_pi(:,:)      =   0.0
    240      nitrcoarse(:,:)      =   0.0
    241      nitrcoarse_pi(:,:)   =   0.0
    242      nitrinscoarse(:,:)   =   0.0
    243      nitrinscoarse_pi(:,:)=   0.0
    244254
    245255  ENDIF !--not aerosol_couple

  • LMDZ5/branches/testing/libf/phylmd/rrtm/suphec.F90

    r2408 r2839  
    129129
    130130IF (LHOOK) CALL DR_HOOK('SUPHEC',0,ZHOOK_HANDLE)
    131 !CALL GSTATS(1811,0) ! MPL 28.11.08
    132 !RVTMP2=RCPV/RCPD-1.0_JPRB   !use cp,moist
    133 RVTMP2=0.0_JPRB              !neglect cp,moist
    134 RHOH2O=RATM/100._JPRB
    135 R2ES=611.21_JPRB*RD/RV
    136 R3LES=17.502_JPRB
    137 R3IES=22.587_JPRB
    138 R4LES=32.19_JPRB
    139 R4IES=-0.7_JPRB
    140 R5LES=R3LES*(RTT-R4LES)
    141 R5IES=R3IES*(RTT-R4IES)
    142 R5ALVCP=R5LES*RLVTT/RCPD
    143 R5ALSCP=R5IES*RLSTT/RCPD
    144 RALVDCP=RLVTT/RCPD
    145 RALSDCP=RLSTT/RCPD
    146 RALFDCP=RLMLT/RCPD
    147 RTWAT=RTT
    148 RTBER=RTT-5._JPRB
    149 RTBERCU=RTT-5.0_JPRB
    150 RTICE=RTT-23._JPRB
    151 RTICECU=RTT-23._JPRB
    152 
    153 RTWAT_RTICE_R=1.0_JPRB/(RTWAT-RTICE)
    154 RTWAT_RTICECU_R=1.0_JPRB/(RTWAT-RTICECU)
    155 IF(NPHYINT == 0) THEN
    156   ISMAX=NSMAX
    157 ELSE
    158   ISMAX=PHYS_GRID%NSMAX
    159 ENDIF
    160 
    161 RKOOP1=2.583_JPRB
    162 RKOOP2=0.48116E-2_JPRB
     131!
     132  IF (OK_BAD_ECMWF_THERMO) THEN
     133!
     134     ! Modify constants defined in suphel.F90 and set RVTMP2 to 0.
     135     ! CALL GSTATS(1811,0) ! MPL 28.11.08
     136     ! RVTMP2=RCPV/RCPD-1.0_JPRB   !use cp,moist
     137     RVTMP2=0.0_JPRB              !neglect cp,moist
     138     RHOH2O=RATM/100._JPRB
     139     R2ES=611.21_JPRB*RD/RV
     140     R3LES=17.502_JPRB
     141     R3IES=22.587_JPRB
     142     R4LES=32.19_JPRB
     143     R4IES=-0.7_JPRB
     144     R5LES=R3LES*(RTT-R4LES)
     145     R5IES=R3IES*(RTT-R4IES)
     146     R5ALVCP=R5LES*RLVTT/RCPD
     147     R5ALSCP=R5IES*RLSTT/RCPD
     148     RALVDCP=RLVTT/RCPD
     149     RALSDCP=RLSTT/RCPD
     150     RALFDCP=RLMLT/RCPD
     151     RTWAT=RTT
     152     RTBER=RTT-5._JPRB
     153     RTBERCU=RTT-5.0_JPRB
     154     RTICE=RTT-23._JPRB
     155     RTICECU=RTT-23._JPRB
     156     
     157     RTWAT_RTICE_R=1.0_JPRB/(RTWAT-RTICE)
     158     RTWAT_RTICECU_R=1.0_JPRB/(RTWAT-RTICECU)
     159     IF(NPHYINT == 0) THEN
     160       ISMAX=NSMAX
     161     ELSE
     162       ISMAX=PHYS_GRID%NSMAX
     163     ENDIF
     164     
     165     RKOOP1=2.583_JPRB
     166     RKOOP2=0.48116E-2_JPRB
     167     
     168  ELSE
     169     ! Keep constants defined in suphel.F90
     170     RTICE=RTT-23._JPRB
     171!
     172  ENDIF  ! (OK_BAD_ECMWF_THERMO)
    163173
    164174!     ------------------------------------------------------------------

  • LMDZ5/branches/testing/libf/phylmd/time_phylmdz_mod.F90

    • Property svn:keywords set to Id
    r2594 r2839  
    11!
    2 ! $Id: $
     2! $Id$
    33!
    44MODULE time_phylmdz_mod
    … …  
    2828    INTEGER,SAVE :: itaufin_phy      ! final iteration (in itau_phy steps)
    2929!$OMP THREADPRIVATE(itaufin_phy)
    30     REAL,SAVE    :: current_time ! current elapsed time (fraction of day) from the begining of the run
     30    REAL,SAVE    :: current_time ! current elapsed time in seconds from the begining of the run
    3131!$OMP THREADPRIVATE(current_time)
    3232   
    … …  
    8484  IMPLICIT NONE
    8585  INCLUDE 'YOMCST.h'
    86     REAL,INTENT(IN) :: pdtphys_
    87     REAL            :: julian_date
    88    
     86  REAL,INTENT(IN) :: pdtphys_
     87  REAL            :: julian_date
     88  INTEGER         :: cur_day
     89  REAL            :: cur_sec
     90
    8991    ! Check if the physics timestep has changed
    9092    IF ( ABS( (pdtphys-pdtphys_) / ((pdtphys+pdtphys_)/2))> 10.*EPSILON(pdtphys_)) THEN
    … …  
    9597   
    9698    ! Update elapsed time since begining of run:
    97     current_time=current_time+pdtphys/rday
     99    current_time = current_time + pdtphys
     100    cur_day = int(current_time/rday)
     101    cur_sec = current_time - (cur_day * rday)
    98102
    99103    ! Compute corresponding Julian date and update calendar
    100     CALL ymds2ju(annee_ref,1,day_ini,(start_time+current_time)*rday,julian_date)
     104    cur_day = cur_day + day_ini
     105    cur_sec = cur_sec + (start_time * rday)
     106    CALL ymds2ju(annee_ref,1, cur_day, cur_sec, julian_date)
    101107    CALL phys_cal_update(julian_date)
    102108   

  • LMDZ5/branches/testing/libf/phylmd/yamada4.F90

    r2729 r2839  
    5757  !                          -> the model can run with longer time-steps.
    5858  !             2016/11/30 (EV etienne.vignon@univ-grenoble-alpes.fr)
    59   !               On met tke (=q2/2) en entrée plutôt que q2
     59  !               On met tke (=q2/2) en entr??e plut??t que q2
    6060  !               On corrige l'update de la tke
    6161  !
    … …  
    120120  DATA first, ipas/.FALSE., 0/
    121121  !$OMP THREADPRIVATE( first,ipas)
    122   LOGICAL hboville
     122  LOGICAL,SAVE :: hboville=.TRUE.
     123  REAL,SAVE :: viscom,viscoh
     124  !$OMP THREADPRIVATE( hboville,viscom,viscoh)
    123125  INTEGER ig, k
    124126  REAL ri, zrif, zalpha, zsm, zsn
    … …  
    130132  REAL zz(klon, klev+1)
    131133  INTEGER iter
     134  REAL dissip(klon,klev), tkeprov,tkeexp, shear(klon,klev), buoy(klon,klev)
    132135  REAL,SAVE :: ric0,ric,rifc, b1, kap
    133136  !$OMP THREADPRIVATE(ric0,ric,rifc,b1,kap)
    … …  
    138141  !$OMP THREADPRIVATE(lmixmin)
    139142  LOGICAL, SAVE :: new_yamada4
    140   !$OMP THREADPRIVATE(new_yamada4)
    141   REAL, SAVE :: yun,ydeux
     143  INTEGER, SAVE :: yamada4_num
     144  !$OMP THREADPRIVATE(new_yamada4,yamada4_num)
     145  REAL, SAVE :: yun,ydeux,disseff
    142146  !$OMP THREADPRIVATE(yun,ydeux)
    143147  REAL frif, falpha, fsm
    … …  
    161165    new_yamada4=.false.
    162166    CALL getin_p('new_yamada4',new_yamada4)
    163 ! Pour garantir la continuite dans la mise au point de CMIP6.
    164 ! Eliminer l'option new_yamada4 des le printemps 2016
     167
    165168    IF (new_yamada4) THEN
     169! Corrections et reglages issus du travail de these d'Etienne Vignon.
    166170       ric=0.143 ! qui donne des valeurs proches des seuils proposes
    167171                 ! dans YAMADA 1983 : sm=0.0845 (0.085 dans Y83)
    … …  
    169173       CALL getin_p('yamada4_ric',ric)
    170174       ric0=0.19489      ! ric=0.195 originalement, mais produisait sm<0
    171        ric=min(ric,ric0) ! Au delà de ric0, sm devient négatif
     175       ric=min(ric,ric0) ! Au dela de ric0, sm devient n??gatif
    172176       rifc=frif(ric)
    173177       seuilsm=fsm(frif(ric))
    … …  
    175179       yun=1.
    176180       ydeux=2.
    177 !      yun=2.
    178 !      ydeux=1.
     181       hboville=.FALSE.
     182       viscom=1.46E-5
     183       viscoh=2.06E-5
     184       lmixmin=0.
     185       yamada4_num=5
    179186    ELSE
    180   !!  DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/
    181187       ric=0.195
    182188       rifc=0.191
    … …  
    185191       yun=2.
    186192       ydeux=1.
     193       hboville=.TRUE.
     194       viscom=0.
     195       viscoh=0.
     196       lmixmin=1.
     197       yamada4_num=0
    187198    ENDIF
     199
    188200    PRINT*,'YAMADA4 RIc, RIfc, Sm_min, Alpha_min',ric,rifc,seuilsm,seuilalpha
    189201    firstcall = .FALSE.
    190     lmixmin=1.
    191202    CALL getin_p('lmixmin',lmixmin)
     203    CALL getin_p('yamada4_hboville',hboville)
     204    CALL getin_p('yamada4_num',yamada4_num)
    192205  END IF
    193206
    … …  
    199212
    200213! On utilise ou non la routine de Holstalg Boville pour les cas tres stables
    201    hboville=.TRUE.
    202214
    203215
    … …  
    290302  !=======================================================================
    291303
    292   ! On commence par calculé q2 à partir de la tke
     304  ! On commence par calculer q2 a partir de la tke
    293305
    294306  IF (new_yamada4) THEN
    … …  
    400412  ELSE IF (iflag_pbl>=10) THEN
    401413
     414    IF (yamada4_num>=1) THEN
     415 
     416    DO k = 2, klev - 1
     417      DO ig=1,ngrid
     418      q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
     419      km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
     420      shear(ig,k)=km(ig, k)*m2(ig, k)
     421      buoy(ig,k)=km(ig, k)*m2(ig, k)*(-1.*rif(ig,k))
     422      dissip(ig,k)=((sqrt(q2(ig,k)))**3)/(b1*l(ig,k))
     423     ENDDO
     424    ENDDO
     425
     426    IF (yamada4_num==1) THEN ! Schema du MAR tel quel
     427       DO k = 2, klev - 1
     428         DO ig=1,ngrid
     429         tkeprov=q2(ig,k)/ydeux
     430         tkeprov= tkeprov*                           &
     431           &  (tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))))/ &
     432           &  (tkeprov+dt*((-1.)*min(0.,buoy(ig,k))+dissip(ig,k)))
     433         q2(ig,k)=tkeprov*ydeux
     434        ENDDO
     435       ENDDO
     436    ELSE IF (yamada4_num==2) THEN ! version modifiee avec integration exacte pour la dissipation
     437       DO k = 2, klev - 1
     438         DO ig=1,ngrid
     439         tkeprov=q2(ig,k)/ydeux
     440         disseff=dissip(ig,k)-min(0.,buoy(ig,k))
     441         tkeprov = tkeprov/(1.+dt*disseff/(2.*tkeprov))**2
     442         tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
     443         q2(ig,k)=tkeprov*ydeux
     444         ! En cas stable, on traite la flotabilite comme la
     445         ! dissipation, en supposant que buoy/q2^3 est constant.
     446         ! Puis on prend la solution exacte
     447        ENDDO
     448       ENDDO
     449    ELSE IF (yamada4_num==3) THEN ! version modifiee avec integration exacte pour la dissipation
     450       DO k = 2, klev - 1
     451         DO ig=1,ngrid
     452         tkeprov=q2(ig,k)/ydeux
     453         disseff=dissip(ig,k)-min(0.,buoy(ig,k))
     454         tkeprov=tkeprov*exp(-dt*disseff/tkeprov)
     455         tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
     456         q2(ig,k)=tkeprov*ydeux
     457         ! En cas stable, on traite la flotabilite comme la
     458         ! dissipation, en supposant que buoy/q2^3 est constant.
     459         ! Puis on prend la solution exacte
     460        ENDDO
     461       ENDDO
     462    ELSE IF (yamada4_num==4) THEN ! version modifiee avec integration exacte pour la dissipation
     463       DO k = 2, klev - 1
     464         DO ig=1,ngrid
     465         tkeprov=q2(ig,k)/ydeux
     466         tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
     467         tkeprov= tkeprov*                           &
     468           &  tkeprov/ &
     469           &  (tkeprov+dt*((-1.)*min(0.,buoy(ig,k))+dissip(ig,k)))
     470         q2(ig,k)=tkeprov*ydeux
     471         ! En cas stable, on traite la flotabilite comme la
     472         ! dissipation, en supposant que buoy/q2^3 est constant.
     473         ! Puis on prend la solution exacte
     474        ENDDO
     475       ENDDO
     476    ELSE IF (yamada4_num==5) THEN ! version modifiee avec integration exacte pour la dissipation
     477       DO k = 2, klev - 1
     478         DO ig=1,ngrid
     479         tkeprov=q2(ig,k)/ydeux
     480         disseff=dissip(ig,k)-min(0.,buoy(ig,k))
     481         tkeexp=exp(-dt*disseff/tkeprov)
     482         tkeprov= shear(ig,k)*tkeprov/disseff*(1.-tkeexp)+tkeprov*tkeexp
     483         q2(ig,k)=tkeprov*ydeux
     484         ! En cas stable, on traite la flotabilite comme la
     485         ! dissipation, en supposant que buoy/q2^3 est constant.
     486         ! Puis on prend la solution exacte
     487        ENDDO
     488       ENDDO
     489    ELSE IF (yamada4_num==6) THEN ! version modifiee avec integration exacte pour la dissipation
     490       DO k = 2, klev - 1
     491         DO ig=1,ngrid
     492         tkeprov=q2(ig,k)/ydeux
     493         tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k),0.)*dt
     494         disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k))
     495         tkeexp=exp(-dt*disseff/tkeprov)
     496         tkeprov= tkeprov*tkeexp
     497         q2(ig,k)=tkeprov*ydeux
     498         ! En cas stable, on traite la flotabilite comme la
     499         ! dissipation, en supposant que buoy/q2^3 est constant.
     500         ! Puis on prend la solution exacte
     501        ENDDO
     502       ENDDO
     503    ENDIF
     504
     505    DO k = 2, klev - 1
     506      DO ig=1,ngrid
     507      q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
     508      ENDDO
     509    ENDDO
     510
     511   ELSE
     512
    402513    DO k = 2, klev - 1
    403514      km(1:ngrid, k) = l(1:ngrid, k)*sqrt(q2(1:ngrid,k))*sm(1:ngrid, k)
    … …  
    410521    END DO
    411522
     523  ENDIF
    412524
    413525  ELSE
    … …  
    529641 END IF ! hboville
    530642
     643! Ajout d'une viscosite moleculaire
     644   km(:,:)=km(:,:)+viscom
     645   kn(:,:)=kn(:,:)+viscoh
     646   kq(:,:)=kq(:,:)+viscoh
     647
    531648  IF (prt_level>1) THEN
    532649    PRINT *, 'YAMADA4 1'
    … …  
    575692
    576693!============================================================================
    577 ! Mise à jour de la tke
     694! Mise a jour de la tke
    578695!============================================================================
    579696
Note: See TracChangeset for help on using the changeset viewer.