Ignore:
Timestamp:
Jul 18, 2016, 9:41:10 PM (8 years ago)
Author:
Laurent Fairhead
Message:

Merged trunk changes r2545:2589 into testing branch

Location:
LMDZ5/branches/testing
Files:
44 edited
3 copied

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing

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

    r2408 r2594  
    8686    INTEGER, ALLOCATABLE, DIMENSION(:) :: hadv  ! index of horizontal trasport schema
    8787    INTEGER, ALLOCATABLE, DIMENSION(:) :: vadv  ! index of vertical trasport schema
     88
     89    INTEGER, ALLOCATABLE, DIMENSION(:) :: hadv_inca  ! index of horizontal trasport schema
     90    INTEGER, ALLOCATABLE, DIMENSION(:) :: vadv_inca  ! index of vertical trasport schema
    8891
    8992    CHARACTER(len=15), ALLOCATABLE, DIMENSION(:) :: tnom_0  ! tracer short name
    … …  
    205208#endif       
    206209       nqtrue=nbtr+nqo
     210
     211       ALLOCATE(hadv_inca(nbtr), vadv_inca(nbtr))
     212
    207213    END IF   ! type_trac
    208214!>jyg
    … …  
    226232!
    227233    ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue),tnom_transp(nqtrue))
     234
    228235!
    229236!jyg<
    … …  
    375382!>jyg
    376383! le module de chimie fournit les noms des traceurs
    377 ! et les schemas d'advection associes.
    378      
     384! et les schemas d'advection associes. excepte pour ceux lus
     385! dans traceur.def
     386       IF (ierr .eq. 0) then
     387          DO iq=1,nqo
     388
     389             write(*,*) 'infotrac 237: iq=',iq
     390             ! CRisi: ajout du nom du fluide transporteur
     391             ! mais rester retro compatible
     392             READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine
     393             write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq)
     394             write(lunout,*) 'tchaine=',trim(tchaine)
     395             write(*,*) 'infotrac 238: IOstatus=',IOstatus
     396             if (IOstatus.ne.0) then
     397                CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1)
     398             endif
     399             ! Y-a-t-il 1 ou 2 noms de traceurs? -> On regarde s'il y a un
     400             ! espace ou pas au milieu de la chaine.
     401             continu=.true.
     402             nouveau_traceurdef=.false.
     403             iiq=1
     404             do while (continu)
     405                if (tchaine(iiq:iiq).eq.' ') then
     406                  nouveau_traceurdef=.true.
     407                  continu=.false.
     408                else if (iiq.lt.LEN_TRIM(tchaine)) then
     409                  iiq=iiq+1
     410                else
     411                  continu=.false.
     412                endif
     413             enddo
     414             write(*,*) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef
     415             if (nouveau_traceurdef) then
     416                write(lunout,*) 'C''est la nouvelle version de traceur.def'
     417                tnom_0(iq)=tchaine(1:iiq-1)
     418                tnom_transp(iq)=tchaine(iiq+1:15)
     419             else
     420                write(lunout,*) 'C''est l''ancienne version de traceur.def'
     421                write(lunout,*) 'On suppose que les traceurs sont tous d''air'
     422                tnom_0(iq)=tchaine
     423                tnom_transp(iq) = 'air'
     424             endif
     425             write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>'
     426             write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>'
     427
     428          END DO !DO iq=1,nqtrue
     429          CLOSE(90) 
     430       ELSE  !! if traceur.def doesn't exist
     431          tnom_0(1)='H2Ov'
     432          tnom_transp(1) = 'air'
     433          tnom_0(2)='H2Ol'
     434          tnom_transp(2) = 'air'
     435          hadv(1) = 10
     436          hadv(2) = 10
     437          vadv(1) = 10
     438          vadv(2) = 10
     439       ENDIF
     440 
    379441#ifdef INCA
    380442       CALL init_transport( &
    381             hadv, &
    382             vadv, &
     443            hadv_inca, &
     444            vadv_inca, &
    383445            conv_flg, &
    384446            pbl_flg,  &
    385447            solsym)
    386448#endif
    387        tnom_0(1)='H2Ov'
    388        tnom_transp(1) = 'air'
    389        tnom_0(2)='H2Ol'
    390        tnom_transp(2) = 'air'
    391        IF (nqo == 3) then
    392           tnom_0(3)='H2Oi'     !! jyg
    393           tnom_transp(3) = 'air'
    394        endif
     449
    395450
    396451!jyg<
    397452       DO iq = nqo+1, nqtrue
     453          hadv(iq) = hadv_inca(iq-nqo)
     454          vadv(iq) = vadv_inca(iq-nqo)
    398455          tnom_0(iq)=solsym(iq-nqo)
    399456          tnom_transp(iq) = 'air'
    400457       END DO
    401 !!       DO iq =3,nqtrue
    402 !!          tnom_0(iq)=solsym(iq-2)
    403 !!       END DO
    404 !!       nqo = 2
    405 !>jyg
    406458
    407459    END IF ! (type_trac == 'inca')

  • LMDZ5/branches/testing/libf/dyn3dmem/dynredem_loc.F90

    r2408 r2594  
    235235
    236236!--- Tracers in file "start_trac.nc" (added by Anne)
     237  lread_inca=.FALSE.
    237238!$OMP MASTER
    238   lread_inca=.FALSE.; fil="start_trac.nc"
     239  fil="start_trac.nc"
    239240  IF(type_trac=='inca') INQUIRE(FILE=fil,EXIST=lread_inca)
    240241  IF(lread_inca) CALL err(NF90_OPEN(fil,NF90_NOWRITE,nid_trac),"open")

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

    r2435 r2594  
    1212                     prad,pg,pr,pcpp,iflag_phys)
    1313  USE dimphy, ONLY: init_dimphy
    14   USE mod_grid_phy_lmdz, ONLY: klon_glo,  & ! number of atmospheric columns (on full grid)
    15                                regular_lonlat  ! regular longitude-latitude grid type
    16   USE mod_phys_lmdz_para, ONLY: klon_omp, & ! number of columns (on local omp grid)
    17                                 klon_omp_begin, & ! start index of local omp subgrid
    18                                 klon_omp_end, & ! end index of local omp subgrid
    19                                 klon_mpi_begin ! start indes of columns (on local mpi grid)
    20   USE geometry_mod, ONLY : init_geometry
     14  USE inigeomphy_mod, ONLY: inigeomphy
     15  USE mod_phys_lmdz_para, ONLY: klon_omp ! number of columns (on local omp grid)
    2116  USE infotrac, ONLY: nqtot, type_trac
    2217  USE infotrac_phy, ONLY: init_infotrac_phy
    23 !  USE comcstphy, ONLY: rradius, & ! planet radius (m)
    24 !                       rr, & ! recuced gas constant: R/molar mass of atm
    25 !                       rg, & ! gravity
    26 !                       rcpp  ! specific heat of the atmosphere
    2718  USE inifis_mod, ONLY: inifis
    2819  USE phyaqua_mod, ONLY: iniaqua
    29   USE physics_distribution_mod, ONLY : init_physics_distribution
    30   USE regular_lonlat_mod, ONLY : init_regular_lonlat, &
    31                                  east, west, north, south, &
    32                                  north_east, north_west, &
    33                                  south_west, south_east
    34   USE mod_interface_dyn_phys, ONLY :  init_interface_dyn_phys
    3520  USE nrtype, ONLY: pi
    3621  IMPLICIT NONE
    … …  
    6954  CHARACTER (LEN=20) :: modname='iniphysiq'
    7055  CHARACTER (LEN=80) :: abort_message
    71   REAL :: total_area_phy, total_area_dyn
    7256
    73   ! boundaries, on global grid
    74   REAL,ALLOCATABLE :: boundslon_reg(:,:)
    75   REAL,ALLOCATABLE :: boundslat_reg(:,:)
    7657
    77   ! global array, on full physics grid:
    78   REAL,ALLOCATABLE :: latfi_glo(:)
    79   REAL,ALLOCATABLE :: lonfi_glo(:)
    80   REAL,ALLOCATABLE :: cufi_glo(:)
    81   REAL,ALLOCATABLE :: cvfi_glo(:)
    82   REAL,ALLOCATABLE :: airefi_glo(:)
    83   REAL,ALLOCATABLE :: boundslonfi_glo(:,:)
    84   REAL,ALLOCATABLE :: boundslatfi_glo(:,:)
     58  ! --> initialize physics distribution, global fields and geometry
     59  ! (i.e. things in phy_common or dynphy_lonlat)
     60  CALL inigeomphy(iim,jjm,nlayer, &
     61               nbp, communicator, &
     62               rlatu,rlatv, &
     63               rlonu,rlonv, &
     64               aire,cu,cv)
    8565
    86   ! local arrays, on given MPI/OpenMP domain:
    87   REAL,ALLOCATABLE,SAVE :: latfi(:)
    88   REAL,ALLOCATABLE,SAVE :: lonfi(:)
    89   REAL,ALLOCATABLE,SAVE :: cufi(:)
    90   REAL,ALLOCATABLE,SAVE :: cvfi(:)
    91   REAL,ALLOCATABLE,SAVE :: airefi(:)
    92   REAL,ALLOCATABLE,SAVE :: boundslonfi(:,:)
    93   REAL,ALLOCATABLE,SAVE :: boundslatfi(:,:)
    94 !$OMP THREADPRIVATE (latfi,lonfi,cufi,cvfi,airefi,boundslonfi,boundslatfi)
    95 
    96   ! Initialize Physics distibution and parameters and interface with dynamics
    97   CALL init_physics_distribution(regular_lonlat,4, &
    98                                  nbp,iim,jjm+1,nlayer,communicator)
    99   CALL init_interface_dyn_phys
    100  
    101   ! init regular global longitude-latitude grid points and boundaries
    102   ALLOCATE(boundslon_reg(iim,2))
    103   ALLOCATE(boundslat_reg(jjm+1,2))
    104  
    105   DO i=1,iim
    106    boundslon_reg(i,east)=rlonu(i)
    107    boundslon_reg(i,west)=rlonu(i+1)
    108   ENDDO
    109 
    110   boundslat_reg(1,north)= PI/2
    111   boundslat_reg(1,south)= rlatv(1)
    112   DO j=2,jjm
    113    boundslat_reg(j,north)=rlatv(j-1)
    114    boundslat_reg(j,south)=rlatv(j)
    115   ENDDO
    116   boundslat_reg(jjm+1,north)= rlatv(jjm)
    117   boundslat_reg(jjm+1,south)= -PI/2
    118 
    119   ! Write values in module regular_lonlat_mod
    120   CALL init_regular_lonlat(iim,jjm+1, rlonv(1:iim), rlatu, &
    121                            boundslon_reg, boundslat_reg)
    122 
    123   ! Generate global arrays on full physics grid
    124   ALLOCATE(latfi_glo(klon_glo),lonfi_glo(klon_glo))
    125   ALLOCATE(cufi_glo(klon_glo),cvfi_glo(klon_glo))
    126   ALLOCATE(airefi_glo(klon_glo))
    127   ALLOCATE(boundslonfi_glo(klon_glo,4))
    128   ALLOCATE(boundslatfi_glo(klon_glo,4))
    129 
    130     ! North pole
    131     latfi_glo(1)=rlatu(1)
    132     lonfi_glo(1)=0.
    133     cufi_glo(1) = cu(1)
    134     cvfi_glo(1) = cv(1)
    135     boundslonfi_glo(1,north_east)=0
    136     boundslatfi_glo(1,north_east)=PI/2
    137     boundslonfi_glo(1,north_west)=2*PI
    138     boundslatfi_glo(1,north_west)=PI/2
    139     boundslonfi_glo(1,south_west)=2*PI
    140     boundslatfi_glo(1,south_west)=rlatv(1)
    141     boundslonfi_glo(1,south_east)=0
    142     boundslatfi_glo(1,south_east)=rlatv(1)
    143     DO j=2,jjm
    144       DO i=1,iim
    145         k=(j-2)*iim+1+i
    146         latfi_glo(k)= rlatu(j)
    147         lonfi_glo(k)= rlonv(i)
    148         cufi_glo(k) = cu((j-1)*(iim+1)+i)
    149         cvfi_glo(k) = cv((j-1)*(iim+1)+i)
    150         boundslonfi_glo(k,north_east)=rlonu(i)
    151         boundslatfi_glo(k,north_east)=rlatv(j-1)
    152         boundslonfi_glo(k,north_west)=rlonu(i+1)
    153         boundslatfi_glo(k,north_west)=rlatv(j-1)
    154         boundslonfi_glo(k,south_west)=rlonu(i+1)
    155         boundslatfi_glo(k,south_west)=rlatv(j)
    156         boundslonfi_glo(k,south_east)=rlonu(i)
    157         boundslatfi_glo(k,south_east)=rlatv(j)
    158       ENDDO
    159     ENDDO
    160     ! South pole
    161     latfi_glo(klon_glo)= rlatu(jjm+1)
    162     lonfi_glo(klon_glo)= 0.
    163     cufi_glo(klon_glo) = cu((iim+1)*jjm+1)
    164     cvfi_glo(klon_glo) = cv((iim+1)*jjm-iim)
    165     boundslonfi_glo(klon_glo,north_east)= 0
    166     boundslatfi_glo(klon_glo,north_east)= rlatv(jjm)
    167     boundslonfi_glo(klon_glo,north_west)= 2*PI
    168     boundslatfi_glo(klon_glo,north_west)= rlatv(jjm)
    169     boundslonfi_glo(klon_glo,south_west)= 2*PI
    170     boundslatfi_glo(klon_glo,south_west)= -PI/2
    171     boundslonfi_glo(klon_glo,south_east)= 0
    172     boundslatfi_glo(klon_glo,south_east)= -Pi/2
    173 
    174     ! build airefi(), mesh area on physics grid
    175     CALL gr_dyn_fi(1,iim+1,jjm+1,klon_glo,aire,airefi_glo)
    176     ! Poles are single points on physics grid
    177     airefi_glo(1)=sum(aire(1:iim,1))
    178     airefi_glo(klon_glo)=sum(aire(1:iim,jjm+1))
    179 
    180     ! Sanity check: do total planet area match between physics and dynamics?
    181     total_area_dyn=sum(aire(1:iim,1:jjm+1))
    182     total_area_phy=sum(airefi_glo(1:klon_glo))
    183     IF (total_area_dyn/=total_area_phy) THEN
    184       WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
    185       WRITE (lunout, *) '     in the dynamics total_area_dyn=', total_area_dyn
    186       WRITE (lunout, *) '  but in the physics total_area_phy=', total_area_phy
    187       IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN
    188         ! stop here if the relative difference is more than 0.001%
    189         abort_message = 'planet total surface discrepancy'
    190         CALL abort_gcm(modname, abort_message, 1)
    191       ENDIF
    192     ENDIF
     66  ! --> now initialize things specific to the phydev physics package
    19367
    19468!$OMP PARALLEL
    195   ! Now generate local lon/lat/cu/cv/area/bounds arrays
    196   ALLOCATE(latfi(klon_omp),lonfi(klon_omp),cufi(klon_omp),cvfi(klon_omp))
    197   ALLOCATE(airefi(klon_omp))
    198   ALLOCATE(boundslonfi(klon_omp,4))
    199   ALLOCATE(boundslatfi(klon_omp,4))
    200 
    201 
    202   offset = klon_mpi_begin - 1
    203   airefi(1:klon_omp) = airefi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    204   cufi(1:klon_omp) = cufi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    205   cvfi(1:klon_omp) = cvfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    206   lonfi(1:klon_omp) = lonfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    207   latfi(1:klon_omp) = latfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    208   boundslonfi(1:klon_omp,:) = boundslonfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
    209   boundslatfi(1:klon_omp,:) = boundslatfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
    210 
    211   ! copy over local grid longitudes and latitudes
    212   CALL init_geometry(klon_omp,lonfi,latfi,boundslonfi,boundslatfi, &
    213                      airefi,cufi,cvfi)
    21469
    21570  ! Initialize physical constants in physics:

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

    r2471 r2594  
    1212                     prad,pg,pr,pcpp,iflag_phys)
    1313  USE dimphy, ONLY: init_dimphy
    14   USE mod_grid_phy_lmdz, ONLY: klon_glo,  & ! number of atmospheric columns (on full grid)
    15                                regular_lonlat, &  ! regular longitude-latitude grid type
    16                                nbp_lon, nbp_lat, nbp_lev
    17   USE mod_phys_lmdz_para, ONLY: klon_omp, & ! number of columns (on local omp grid)
    18                                 klon_omp_begin, & ! start index of local omp subgrid
    19                                 klon_omp_end, & ! end index of local omp subgrid
    20                                 klon_mpi_begin ! start indes of columns (on local mpi grid)
    21   USE geometry_mod, ONLY : init_geometry
     14  USE inigeomphy_mod, ONLY: inigeomphy
     15  USE mod_grid_phy_lmdz, ONLY: klon_glo ! number of atmospheric columns (on full grid)
     16  USE mod_phys_lmdz_para, ONLY: klon_omp ! number of columns (on local omp grid)
    2217  USE vertical_layers_mod, ONLY : init_vertical_layers
    2318  USE infotrac, ONLY: nqtot,nqo,nbtr,tname,ttext,type_trac,&
    … …  
    3934  USE phystokenc_mod, ONLY: init_phystokenc
    4035  USE phyaqua_mod, ONLY: iniaqua
    41   USE physics_distribution_mod, ONLY : init_physics_distribution
    42   USE regular_lonlat_mod, ONLY : init_regular_lonlat, &
    43                                  east, west, north, south, &
    44                                  north_east, north_west, &
    45                                  south_west, south_east
    46   USE mod_interface_dyn_phys, ONLY :  init_interface_dyn_phys
    4736#ifdef INCA
    4837  USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic
    … …  
    9281  CHARACTER (LEN=20) :: modname = 'iniphysiq'
    9382  CHARACTER (LEN=80) :: abort_message
    94   REAL :: total_area_phy, total_area_dyn
    9583
    96   ! boundaries, on global grid
    97   REAL,ALLOCATABLE :: boundslon_reg(:,:)
    98   REAL,ALLOCATABLE :: boundslat_reg(:,:)
    99 
    100   ! global array, on full physics grid:
    101   REAL,ALLOCATABLE :: latfi_glo(:)
    102   REAL,ALLOCATABLE :: lonfi_glo(:)
    103   REAL,ALLOCATABLE :: cufi_glo(:)
    104   REAL,ALLOCATABLE :: cvfi_glo(:)
    105   REAL,ALLOCATABLE :: airefi_glo(:)
    106   REAL,ALLOCATABLE :: boundslonfi_glo(:,:)
    107   REAL,ALLOCATABLE :: boundslatfi_glo(:,:)
    108 
    109   ! local arrays, on given MPI/OpenMP domain:
    110   REAL,ALLOCATABLE,SAVE :: latfi(:)
    111   REAL,ALLOCATABLE,SAVE :: lonfi(:)
    112   REAL,ALLOCATABLE,SAVE :: cufi(:)
    113   REAL,ALLOCATABLE,SAVE :: cvfi(:)
    114   REAL,ALLOCATABLE,SAVE :: airefi(:)
    115   REAL,ALLOCATABLE,SAVE :: boundslonfi(:,:)
    116   REAL,ALLOCATABLE,SAVE :: boundslatfi(:,:)
    117 !$OMP THREADPRIVATE (latfi,lonfi,cufi,cvfi,airefi,boundslonfi,boundslatfi)
    11884
    11985#ifndef CPP_PARA
    … …  
    12389#endif
    12490
    125   ! Initialize Physics distibution and parameters and interface with dynamics
    126   IF (ii*jj>1) THEN ! general 3D case
    127     CALL init_physics_distribution(regular_lonlat,4, &
    128                                  nbp,ii,jj+1,nlayer,communicator)
    129   ELSE ! For 1D model
    130     CALL init_physics_distribution(regular_lonlat,4, &
    131                                  1,1,1,nlayer,communicator) 
    132   ENDIF
    133   CALL init_interface_dyn_phys
     91  ! --> initialize physics distribution, global fields and geometry
     92  ! (i.e. things in phy_common or dynphy_lonlat)
     93  CALL inigeomphy(ii,jj,nlayer, &
     94               nbp, communicator, &
     95               rlatu,rlatv, &
     96               rlonu,rlonv, &
     97               aire,cu,cv)
     98
     99  ! --> now initialize things specific to the phylmd physics package
    134100 
    135   ! init regular global longitude-latitude grid points and boundaries
    136   ALLOCATE(boundslon_reg(ii,2))
    137   ALLOCATE(boundslat_reg(jj+1,2))
    138  
    139   DO i=1,ii
    140    boundslon_reg(i,east)=rlonu(i)
    141    boundslon_reg(i,west)=rlonu(i+1)
    142   ENDDO
    143 
    144   boundslat_reg(1,north)= PI/2
    145   boundslat_reg(1,south)= rlatv(1)
    146   DO j=2,jj
    147    boundslat_reg(j,north)=rlatv(j-1)
    148    boundslat_reg(j,south)=rlatv(j)
    149   ENDDO
    150   boundslat_reg(jj+1,north)= rlatv(jj)
    151   boundslat_reg(jj+1,south)= -PI/2
    152 
    153   ! Write values in module regular_lonlat_mod
    154   CALL init_regular_lonlat(ii,jj+1, rlonv(1:ii), rlatu, &
    155                            boundslon_reg, boundslat_reg)
    156 
    157   ! Generate global arrays on full physics grid
    158   ALLOCATE(latfi_glo(klon_glo),lonfi_glo(klon_glo))
    159   ALLOCATE(cufi_glo(klon_glo),cvfi_glo(klon_glo))
    160   ALLOCATE(airefi_glo(klon_glo))
    161   ALLOCATE(boundslonfi_glo(klon_glo,4))
    162   ALLOCATE(boundslatfi_glo(klon_glo,4))
    163  
    164 
    165   IF (klon_glo>1) THEN ! general case
    166     ! North pole
    167     latfi_glo(1)=rlatu(1)
    168     lonfi_glo(1)=0.
    169     cufi_glo(1) = cu(1)
    170     cvfi_glo(1) = cv(1)
    171     boundslonfi_glo(1,north_east)=0
    172     boundslatfi_glo(1,north_east)=PI/2
    173     boundslonfi_glo(1,north_west)=2*PI
    174     boundslatfi_glo(1,north_west)=PI/2
    175     boundslonfi_glo(1,south_west)=2*PI
    176     boundslatfi_glo(1,south_west)=rlatv(1)
    177     boundslonfi_glo(1,south_east)=0
    178     boundslatfi_glo(1,south_east)=rlatv(1)
    179     DO j=2,jj
    180       DO i=1,ii
    181         k=(j-2)*ii+1+i
    182         latfi_glo(k)= rlatu(j)
    183         lonfi_glo(k)= rlonv(i)
    184         cufi_glo(k) = cu((j-1)*(ii+1)+i)
    185         cvfi_glo(k) = cv((j-1)*(ii+1)+i)
    186         boundslonfi_glo(k,north_east)=rlonu(i)
    187         boundslatfi_glo(k,north_east)=rlatv(j-1)
    188         boundslonfi_glo(k,north_west)=rlonu(i+1)
    189         boundslatfi_glo(k,north_west)=rlatv(j-1)
    190         boundslonfi_glo(k,south_west)=rlonu(i+1)
    191         boundslatfi_glo(k,south_west)=rlatv(j)
    192         boundslonfi_glo(k,south_east)=rlonu(i)
    193         boundslatfi_glo(k,south_east)=rlatv(j)
    194       ENDDO
    195     ENDDO
    196     ! South pole
    197     latfi_glo(klon_glo)= rlatu(jj+1)
    198     lonfi_glo(klon_glo)= 0.
    199     cufi_glo(klon_glo) = cu((ii+1)*jj+1)
    200     cvfi_glo(klon_glo) = cv((ii+1)*jj-ii)
    201     boundslonfi_glo(klon_glo,north_east)= 0
    202     boundslatfi_glo(klon_glo,north_east)= rlatv(jj)
    203     boundslonfi_glo(klon_glo,north_west)= 2*PI
    204     boundslatfi_glo(klon_glo,north_west)= rlatv(jj)
    205     boundslonfi_glo(klon_glo,south_west)= 2*PI
    206     boundslatfi_glo(klon_glo,south_west)= -PI/2
    207     boundslonfi_glo(klon_glo,south_east)= 0
    208     boundslatfi_glo(klon_glo,south_east)= -Pi/2
    209 
    210     ! build airefi_glo(), mesh area on physics grid
    211     CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,aire,airefi_glo)
    212     ! Poles are single points on physics grid
    213     airefi_glo(1)=sum(aire(1:ii,1))
    214     airefi_glo(klon_glo)=sum(aire(1:ii,jj+1))
    215 
    216     ! Sanity check: do total planet area match between physics and dynamics?
    217     total_area_dyn=sum(aire(1:ii,1:jj+1))
    218     total_area_phy=sum(airefi_glo(1:klon_glo))
    219     IF (total_area_dyn/=total_area_phy) THEN
    220       WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
    221       WRITE (lunout, *) '     in the dynamics total_area_dyn=', total_area_dyn
    222       WRITE (lunout, *) '  but in the physics total_area_phy=', total_area_phy
    223       IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN
    224         ! stop here if the relative difference is more than 0.001%
    225         abort_message = 'planet total surface discrepancy'
    226         CALL abort_gcm(modname, abort_message, 1)
    227       ENDIF
    228     ENDIF
    229   ELSE ! klon_glo==1, running the 1D model
    230     ! just copy over input values
    231     latfi_glo(1)=rlatu(1)
    232     lonfi_glo(1)=rlonv(1)
    233     cufi_glo(1)=cu(1)
    234     cvfi_glo(1)=cv(1)
    235     airefi_glo(1)=aire(1,1)
    236     boundslonfi_glo(1,north_east)=rlonu(1)
    237     boundslatfi_glo(1,north_east)=PI/2
    238     boundslonfi_glo(1,north_west)=rlonu(2)
    239     boundslatfi_glo(1,north_west)=PI/2
    240     boundslonfi_glo(1,south_west)=rlonu(2)
    241     boundslatfi_glo(1,south_west)=rlatv(1)
    242     boundslonfi_glo(1,south_east)=rlonu(1)
    243     boundslatfi_glo(1,south_east)=rlatv(1)
    244   ENDIF ! of IF (klon_glo>1)
    245 
    246101!$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/)
    247   ! Now generate local lon/lat/cu/cv/area/bounds arrays
    248   ALLOCATE(latfi(klon_omp),lonfi(klon_omp),cufi(klon_omp),cvfi(klon_omp))
    249   ALLOCATE(airefi(klon_omp))
    250   ALLOCATE(boundslonfi(klon_omp,4))
    251   ALLOCATE(boundslatfi(klon_omp,4))
    252 
    253 
    254   offset = klon_mpi_begin - 1
    255   airefi(1:klon_omp) = airefi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    256   cufi(1:klon_omp) = cufi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    257   cvfi(1:klon_omp) = cvfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    258   lonfi(1:klon_omp) = lonfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    259   latfi(1:klon_omp) = latfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    260   boundslonfi(1:klon_omp,:) = boundslonfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
    261   boundslatfi(1:klon_omp,:) = boundslatfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
    262 
    263   ! copy over local grid longitudes and latitudes
    264   CALL init_geometry(klon_omp,lonfi,latfi,boundslonfi,boundslatfi, &
    265                      airefi,cufi,cvfi)
    266102
    267103  ! copy over preff , ap(), bp(), etc

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

    r2542 r2594  
    117117    ALLOCATE(pctsrf(klon,nbsrf))
    118118    CALL start_init_subsurf(.FALSE.)
     119  !--- TO MATCH EXACTLY WHAT WOULD BE DONE IN etat0phys_netcdf
     120    WHERE(   masque(:,:)<EPSFRA) masque(:,:)=0.
     121    WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
    119122  END IF
    120123
    … …  
    336339  REAL, ALLOCATABLE :: champan(:,:,:)
    337340!--- input files
     341  CHARACTER(LEN=20) :: fnam_m, fnam_p     ! previous/next files names
    338342  CHARACTER(LEN=20) :: cal_in             ! calendar
    339343  CHARACTER(LEN=20) :: unit_sic           ! attribute unit in sea-ice file
    … …  
    345349  LOGICAL           :: extrp              ! flag for extrapolation
    346350  REAL              :: chmin, chmax
    347   INTEGER ierr
     351  INTEGER ierr, idx
    348352  integer n_extrap ! number of extrapolated points
    349353  logical skip
    … …  
    360364  END SELECT
    361365  extrp=.FALSE.; IF(PRESENT(flag).AND.mode=='SST') extrp=flag
     366  idx=INDEX(fnam,'.nc')-1
    362367
    363368!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
    … …  
    393398!--- Time (variable is not needed - it is rebuilt - but calendar is)
    394399  CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep), fnam)
    395   ALLOCATE(timeyear(MAX(lmdep,14)))
     400  ALLOCATE(timeyear(lmdep+2))
    396401  CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
    397402  cal_in=' '
    … …  
    412417  IF(lmdep==12) THEN
    413418    timeyear=mid_month(anneeref, cal_in, ndays_in)
    414     CALL msg(0,'Monthly periodic input file (perpetual run).')
    415   ELSE IF(lmdep==14) THEN
    416     timeyear=mid_month(anneeref, cal_in, ndays_in)
    417     CALL msg(0,'Monthly 14-records input file (interannual run).')
     419    CALL msg(0,'Monthly input file(s) for '//TRIM(title)//'.')
    418420  ELSE IF(lmdep==ndays_in) THEN
    419     timeyear=[(REAL(k)-0.5,k=1,ndays_in)]
     421    timeyear=[(REAL(k)-0.5,k=0,ndays_in+1)]
    420422    CALL msg(0,'Daily input file (no time interpolation).')
    421423  ELSE
    422424    WRITE(mess,'(a,i3,a,i3,a)')'Mismatching input file: found',lmdep,        &
    423       ' records, 12/14/',ndays_in,' (periodic/interannual/daily) needed'
     425      ' records, 12/',ndays_in,' (monthly/daily needed).'
    424426    CALL abort_physic('mid_months',TRIM(mess),1)
    425427  END IF
    426428
    427429!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
    428   ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, MAX(lmdep,14)))
     430  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep+2))
    429431  IF(extrp) ALLOCATE(work(imdep, jmdep))
    430432  CALL msg(5,'')
    … …  
    446448      WHERE(NINT(mask)/=1) champint=0.001
    447449    END IF
    448     !--- Special case for periodic input file: index shifted
    449     ll = l; IF(lmdep==12) ll = l + 1
    450     champtime(:, :, ll)=champint
     450    champtime(:, :, l+1)=champint
    451451  END DO
    452   IF(lmdep==12) THEN
    453     champtime(:,:, 1)=champtime(:,:,13)
    454     champtime(:,:,14)=champtime(:,:, 2)
    455   END IF
    456452  CALL ncerr(NF90_CLOSE(ncid), fnam)
    457453
     454!--- FIRST RECORD: LAST ONE OF PREVIOUS YEAR (CURRENT YEAR IF UNAVAILABLE)
     455  fnam_m=fnam(1:idx)//'_m.nc'
     456  IF(NF90_OPEN(fnam_m,NF90_NOWRITE,ncid)==NF90_NOERR) THEN
     457    CALL msg(0,'Reading previous year file ("'//TRIM(fnam_m)//'") last record for '//TRIM(title))
     458    CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_m)
     459    CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam_m)
     460    CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), len=l), fnam_m)
     461    CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam_m)
     462    CALL ncerr(NF90_CLOSE(ncid), fnam_m)
     463    CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
     464    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
     465    IF(mode=='RUG') champ=LOG(champ)
     466    CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
     467    IF(mode=='RUG') THEN
     468      champint=EXP(champint)
     469      WHERE(NINT(mask)/=1) champint=0.001
     470    END IF
     471    champtime(:, :, 1)=champint
     472  ELSE
     473    CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") last record for '//TRIM(title))
     474    champtime(:, :, 1)=champtime(:, :, lmdep+1)
     475  END IF
     476
     477!--- LAST RECORD: FIRST ONE OF NEXT YEAR (CURRENT YEAR IF UNAVAILABLE)
     478  fnam_p=fnam(1:idx)//'_p.nc'
     479  IF(NF90_OPEN(fnam_p,NF90_NOWRITE,ncid)==NF90_NOERR) THEN
     480    CALL msg(0,'Reading previous year file ("'//TRIM(fnam_p)//'") first record for '//TRIM(title))
     481    CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_p)
     482    CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,1],[imdep,jmdep,1]),fnam_p)
     483    CALL ncerr(NF90_CLOSE(ncid), fnam_p)
     484    CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
     485    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
     486    IF(mode=='RUG') champ=LOG(champ)
     487    CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
     488    IF(mode=='RUG') THEN
     489      champint=EXP(champint)
     490      WHERE(NINT(mask)/=1) champint=0.001
     491    END IF
     492    champtime(:, :, lmdep+2)=champint
     493  ELSE
     494    CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") first record for '//TRIM(title))
     495    champtime(:, :, lmdep+2)=champtime(:, :, 2)
     496  END IF
    458497  DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
    459498  IF(extrp) DEALLOCATE(work)
    … …  
    477516     END IF
    478517  END IF
    479   ALLOCATE(yder(MAX(lmdep,14)), champan(iip1, jjp1, ndays))
     518  ALLOCATE(yder(lmdep+2), champan(iip1, jjp1, ndays))
    480519  IF(lmdep==ndays_in) THEN
    481520     champan(1:iim,:,:)=champtime
    … …  
    546585!
    547586!-------------------------------------------------------------------------------
     587  USE grid_noro_m, ONLY: grid_noro0
    548588  IMPLICIT NONE
    549589!===============================================================================
    … …  
    599639
    600640END SUBROUTINE start_init_orog0
    601 !
    602 !-------------------------------------------------------------------------------
    603 
    604 
    605 !-------------------------------------------------------------------------------
    606 !
    607 SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)
    608 !
    609 !===============================================================================
    610 ! Purpose: Extracted from grid_noro to provide geopotential height for dynamics
    611 !          without any call to physics subroutines.
    612 !===============================================================================
    613   IMPLICIT NONE
    614 !-------------------------------------------------------------------------------
    615 ! Arguments:
    616   REAL, INTENT(IN)   :: xd(:), yd(:) !--- INPUT  COORDINATES     (imdp) (jmdp)
    617   REAL, INTENT(IN)   :: zd(:,:)      !--- INPUT  FIELD           (imdp,jmdp)
    618   REAL, INTENT(IN)   :: x(:), y(:)   !--- OUTPUT COORDINATES     (imar+1) (jmar)
    619   REAL, INTENT(OUT)  :: zphi(:,:)    !--- GEOPOTENTIAL           (imar+1,jmar)
    620   REAL, INTENT(INOUT):: mask(:,:)    !--- MASK                   (imar+1,jmar)
    621 !-------------------------------------------------------------------------------
    622 ! Local variables:
    623   CHARACTER(LEN=256) :: modname="grid_noro0"
    624   REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
    625   REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
    626   REAL, ALLOCATABLE :: weight(:,:)                ! dim (imar+1,jmar)
    627   REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:)   ! dim (imar+1,jmar)
    628   REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
    629   REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
    630   REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
    631   LOGICAL :: masque_lu
    632   INTEGER :: i, ii, imdp, imar, iext
    633   INTEGER :: j, jj, jmdp, jmar, nn
    634   REAL    :: xpi, zlenx, weighx, xincr,  zbordnor, zmeanor, zweinor, zbordest
    635   REAL    :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue
    636 !-------------------------------------------------------------------------------
    637   imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
    638   jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
    639   imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
    640   jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
    641   IF(imar/=iim)   CALL abort_gcm(TRIM(modname),'imar/=iim'  ,1)
    642   IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
    643   iext=imdp/10
    644   xpi = ACOS(-1.)
    645   rad = 6371229.
    646 
    647 !--- ARE WE USING A READ MASK ?
    648   masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
    649   WRITE(lunout,*)'Masque lu: ',masque_lu
    650 
    651 !--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
    652   ALLOCATE(xusn(imdp+2*iext))
    653   xusn(1     +iext:imdp  +iext)=xd(:)
    654   xusn(1          :       iext)=xd(1+imdp-iext:imdp)-2.*xpi
    655   xusn(1+imdp+iext:imdp+2*iext)=xd(1          :iext)+2.*xpi
    656 
    657   ALLOCATE(yusn(jmdp+2))
    658   yusn(1       )=yd(1)   +(yd(1)   -yd(2))
    659   yusn(2:jmdp+1)=yd(:)
    660   yusn(  jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
    661 
    662   ALLOCATE(zusn(imdp+2*iext,jmdp+2))
    663   zusn(1       +iext:imdp  +iext,2:jmdp+1)=zd  (:                   ,     :)
    664   zusn(1            :       iext,2:jmdp+1)=zd  (imdp-iext+1:imdp    ,     :)
    665   zusn(1+imdp  +iext:imdp+2*iext,2:jmdp+1)=zd  (1:iext              ,     :)
    666   zusn(1            :imdp/2+iext,       1)=zusn(1+imdp/2:imdp  +iext,     2)
    667   zusn(1+imdp/2+iext:imdp+2*iext,       1)=zusn(1       :imdp/2+iext,     2)
    668   zusn(1            :imdp/2+iext,  jmdp+2)=zusn(1+imdp/2:imdp  +iext,jmdp+1)
    669   zusn(1+imdp/2+iext:imdp+2*iext,  jmdp+2)=zusn(1       :imdp/2+iext,jmdp+1)
    670 
    671 !--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
    672   ALLOCATE(a(imar+1),b(imar+1))
    673   b(1:imar)=(x(1:imar  )+ x(2:imar+1))/2.0
    674   b(imar+1)= x(  imar+1)+(x(  imar+1)-x(imar))/2.0
    675   a(1)=x(1)-(x(2)-x(1))/2.0
    676   a(2:imar+1)= b(1:imar)
    677 
    678   ALLOCATE(c(jmar),d(jmar))
    679   d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
    680   d(  jmar  )= y(  jmar  )+(y(  jmar)-y(jmar-1))/2.0
    681   c(1)=y(1)-(y(2)-y(1))/2.0
    682   c(2:jmar)=d(1:jmar-1)
    683 
    684 !--- INITIALIZATIONS:
    685   ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
    686   ALLOCATE(zmea  (imar+1,jmar)); zmea  (:,:)= 0.0
    687 
    688 !--- SUMMATION OVER GRIDPOINT AREA
    689   zleny=xpi/REAL(jmdp)*rad
    690   xincr=xpi/REAL(jmdp)/2.
    691   ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
    692   ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
    693   DO ii = 1, imar+1
    694     DO jj = 1, jmar
    695       DO j = 2,jmdp+1
    696         zlenx  =zleny  *COS(yusn(j))
    697         zbordnor=(xincr+c(jj)-yusn(j))*rad
    698         zbordsud=(xincr-d(jj)+yusn(j))*rad
    699         weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
    700         IF(weighy/=0) THEN
    701           DO i = 2, imdp+2*iext-1
    702             zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))
    703             zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))
    704             weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
    705             IF(weighx/=0)THEN
    706               num_tot(ii,jj)=num_tot(ii,jj)+1.0
    707               IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
    708               weight(ii,jj)=weight(ii,jj)+weighx*weighy
    709               zmea  (ii,jj)=zmea  (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
    710             END IF
    711           END DO
    712         END IF
    713       END DO
    714     END DO
    715   END DO
    716 
    717 !--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
    718   IF(.NOT.masque_lu) THEN
    719     WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
    720   END IF
    721   nn=COUNT(weight(:,1:jmar-1)==0.0)
    722   IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
    723   WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
    724 
    725 !--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)
    726   ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
    727   WHERE(mask>=0.1) mask_tmp = 1.
    728   WHERE(weight(:,:)/=0.0)
    729     zphi(:,:)=mask_tmp(:,:)*zmea(:,:)
    730     zmea(:,:)=mask_tmp(:,:)*zmea(:,:)
    731   END WHERE
    732   WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
    733 
    734 !--- Values at poles
    735   zphi(imar+1,:)=zphi(1,:)
    736 
    737   zweinor=SUM(weight(1:imar,   1),DIM=1)
    738   zweisud=SUM(weight(1:imar,jmar),DIM=1)
    739   zmeanor=SUM(weight(1:imar,   1)*zmea(1:imar,   1),DIM=1)
    740   zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)
    741   zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud
    742 
    743 END SUBROUTINE grid_noro0
    744641!
    745642!-------------------------------------------------------------------------------

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

    r2435 r2594  
    1212                     prad,pg,pr,pcpp,iflag_phys)
    1313  USE dimphy, ONLY: init_dimphy
    14   USE mod_grid_phy_lmdz, ONLY: klon_glo,  & ! number of atmospheric columns (on full grid)
    15                                regular_lonlat  ! regular longitude-latitude grid type
    16   USE mod_phys_lmdz_para, ONLY: klon_omp, & ! number of columns (on local omp grid)
    17                                 klon_omp_begin, & ! start index of local omp subgrid
    18                                 klon_omp_end, & ! end index of local omp subgrid
    19                                 klon_mpi_begin ! start indes of columns (on local mpi grid)
    20   USE geometry_mod, ONLY : init_geometry
     14  USE inigeomphy_mod, ONLY: inigeomphy
    2115  USE infotrac, ONLY: nqtot
    2216  USE comcstphy, ONLY: rradius, & ! planet radius (m)
    … …  
    2418                       rg, & ! gravity
    2519                       rcpp  ! specific heat of the atmosphere
    26 !  USE phyaqua_mod, ONLY: iniaqua
    27   USE physics_distribution_mod, ONLY : init_physics_distribution
    28   USE regular_lonlat_mod, ONLY : init_regular_lonlat, &
    29                                  east, west, north, south, &
    30                                  north_east, north_west, &
    31                                  south_west, south_east
    32   USE mod_interface_dyn_phys, ONLY :  init_interface_dyn_phys
    3320  USE infotrac_phy, ONLY: init_infotrac_phy
    3421  USE nrtype, ONLY: pi
    … …  
    6855  CHARACTER (LEN=20) :: modname='iniphysiq'
    6956  CHARACTER (LEN=80) :: abort_message
    70   REAL :: total_area_phy, total_area_dyn
    71 
    72   ! boundaries, on global grid
    73   REAL,ALLOCATABLE :: boundslon_reg(:,:)
    74   REAL,ALLOCATABLE :: boundslat_reg(:,:)
    75 
    76   ! global array, on full physics grid:
    77   REAL,ALLOCATABLE :: latfi_glo(:)
    78   REAL,ALLOCATABLE :: lonfi_glo(:)
    79   REAL,ALLOCATABLE :: cufi_glo(:)
    80   REAL,ALLOCATABLE :: cvfi_glo(:)
    81   REAL,ALLOCATABLE :: airefi_glo(:)
    82   REAL,ALLOCATABLE :: boundslonfi_glo(:,:)
    83   REAL,ALLOCATABLE :: boundslatfi_glo(:,:)
    84 
    85   ! local arrays, on given MPI/OpenMP domain:
    86   REAL,ALLOCATABLE,SAVE :: latfi(:)
    87   REAL,ALLOCATABLE,SAVE :: lonfi(:)
    88   REAL,ALLOCATABLE,SAVE :: cufi(:)
    89   REAL,ALLOCATABLE,SAVE :: cvfi(:)
    90   REAL,ALLOCATABLE,SAVE :: airefi(:)
    91   REAL,ALLOCATABLE,SAVE :: boundslonfi(:,:)
    92   REAL,ALLOCATABLE,SAVE :: boundslatfi(:,:)
    93 !$OMP THREADPRIVATE (latfi,lonfi,cufi,cvfi,airefi,boundslonfi,boundslatfi)
    94 
    95   ! Initialize Physics distibution and parameters and interface with dynamics
    96   CALL init_physics_distribution(regular_lonlat,4, &
    97                                  nbp,ii,jj+1,nlayer,communicator)
    98   CALL init_interface_dyn_phys
    99  
    100   ! init regular global longitude-latitude grid points and boundaries
    101   ALLOCATE(boundslon_reg(ii,2))
    102   ALLOCATE(boundslat_reg(jj+1,2))
    103  
    104   DO i=1,ii
    105    boundslon_reg(i,east)=rlonu(i)
    106    boundslon_reg(i,west)=rlonu(i+1)
    107   ENDDO
    108 
    109   boundslat_reg(1,north)= PI/2
    110   boundslat_reg(1,south)= rlatv(1)
    111   DO j=2,jj
    112    boundslat_reg(j,north)=rlatv(j-1)
    113    boundslat_reg(j,south)=rlatv(j)
    114   ENDDO
    115   boundslat_reg(jj+1,north)= rlatv(jj)
    116   boundslat_reg(jj+1,south)= -PI/2
    117 
    118   ! Write values in module regular_lonlat_mod
    119   CALL init_regular_lonlat(ii,jj+1, rlonv(1:ii), rlatu, &
    120                            boundslon_reg, boundslat_reg)
    121 
    122   ! Generate global arrays on full physics grid
    123   ALLOCATE(latfi_glo(klon_glo),lonfi_glo(klon_glo))
    124   ALLOCATE(cufi_glo(klon_glo),cvfi_glo(klon_glo))
    125   ALLOCATE(airefi_glo(klon_glo))
    126   ALLOCATE(boundslonfi_glo(klon_glo,4))
    127   ALLOCATE(boundslatfi_glo(klon_glo,4))
    128 
    129   IF (klon_glo>1) THEN ! general case
    130     ! North pole
    131     latfi_glo(1)=rlatu(1)
    132     lonfi_glo(1)=0.
    133     cufi_glo(1) = cu(1)
    134     cvfi_glo(1) = cv(1)
    135     boundslonfi_glo(1,north_east)=0
    136     boundslatfi_glo(1,north_east)=PI/2
    137     boundslonfi_glo(1,north_west)=2*PI
    138     boundslatfi_glo(1,north_west)=PI/2
    139     boundslonfi_glo(1,south_west)=2*PI
    140     boundslatfi_glo(1,south_west)=rlatv(1)
    141     boundslonfi_glo(1,south_east)=0
    142     boundslatfi_glo(1,south_east)=rlatv(1)
    143     DO j=2,jj
    144       DO i=1,ii
    145         k=(j-2)*ii+1+i
    146         latfi_glo((j-2)*ii+1+i)= rlatu(j)
    147         lonfi_glo((j-2)*ii+1+i)= rlonv(i)
    148         cufi_glo((j-2)*ii+1+i) = cu((j-1)*(ii+1)+i)
    149         cvfi_glo((j-2)*ii+1+i) = cv((j-1)*(ii+1)+i)
    150         boundslonfi_glo(k,north_east)=rlonu(i)
    151         boundslatfi_glo(k,north_east)=rlatv(j-1)
    152         boundslonfi_glo(k,north_west)=rlonu(i+1)
    153         boundslatfi_glo(k,north_west)=rlatv(j-1)
    154         boundslonfi_glo(k,south_west)=rlonu(i+1)
    155         boundslatfi_glo(k,south_west)=rlatv(j)
    156         boundslonfi_glo(k,south_east)=rlonu(i)
    157         boundslatfi_glo(k,south_east)=rlatv(j)
    158       ENDDO
    159     ENDDO
    160     ! South pole
    161     latfi_glo(klon_glo)= rlatu(jj+1)
    162     lonfi_glo(klon_glo)= 0.
    163     cufi_glo(klon_glo) = cu((ii+1)*jj+1)
    164     cvfi_glo(klon_glo) = cv((ii+1)*jj-ii)
    165     boundslonfi_glo(klon_glo,north_east)= 0
    166     boundslatfi_glo(klon_glo,north_east)= rlatv(jj)
    167     boundslonfi_glo(klon_glo,north_west)= 2*PI
    168     boundslatfi_glo(klon_glo,north_west)= rlatv(jj)
    169     boundslonfi_glo(klon_glo,south_west)= 2*PI
    170     boundslatfi_glo(klon_glo,south_west)= -PI/2
    171     boundslonfi_glo(klon_glo,south_east)= 0
    172     boundslatfi_glo(klon_glo,south_east)= -Pi/2
    173 
    174     ! build airefi(), mesh area on physics grid
    175     CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,aire,airefi_glo)
    176     ! Poles are single points on physics grid
    177     airefi_glo(1)=sum(aire(1:ii,1))
    178     airefi_glo(klon_glo)=sum(aire(1:ii,jj+1))
    179 
    180     ! Sanity check: do total planet area match between physics and dynamics?
    181     total_area_dyn=sum(aire(1:ii,1:jj+1))
    182     total_area_phy=sum(airefi_glo(1:klon_glo))
    183     IF (total_area_dyn/=total_area_phy) THEN
    184       WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
    185       WRITE (lunout, *) '     in the dynamics total_area_dyn=', total_area_dyn
    186       WRITE (lunout, *) '  but in the physics total_area_phy=', total_area_phy
    187       IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN
    188         ! stop here if the relative difference is more than 0.001%
    189         abort_message = 'planet total surface discrepancy'
    190         CALL abort_gcm(modname, abort_message, 1)
    191       ENDIF
    192     ENDIF
    193   ELSE ! klon_glo==1, running the 1D model
    194     ! just copy over input values
    195     latfi_glo(1)=rlatu(1)
    196     lonfi_glo(1)=rlonv(1)
    197     cufi_glo(1)=cu(1)
    198     cvfi_glo(1)=cv(1)
    199     airefi_glo(1)=aire(1,1)
    200     boundslonfi_glo(1,north_east)=rlonu(1)
    201     boundslatfi_glo(1,north_east)=PI/2
    202     boundslonfi_glo(1,north_west)=rlonu(2)
    203     boundslatfi_glo(1,north_west)=PI/2
    204     boundslonfi_glo(1,south_west)=rlonu(2)
    205     boundslatfi_glo(1,south_west)=rlatv(1)
    206     boundslonfi_glo(1,south_east)=rlonu(1)
    207     boundslatfi_glo(1,south_east)=rlatv(1)
    208   ENDIF ! of IF (klon_glo>1)
    209 
    210 !$OMP PARALLEL
    211   ! Now generate local lon/lat/cu/cv/area/bounds arrays
    212   ALLOCATE(latfi(klon_omp),lonfi(klon_omp),cufi(klon_omp),cvfi(klon_omp))
    213   ALLOCATE(airefi(klon_omp))
    214   ALLOCATE(boundslonfi(klon_omp,4))
    215   ALLOCATE(boundslatfi(klon_omp,4))
    21657
    21758
    218   offset = klon_mpi_begin - 1
    219   airefi(1:klon_omp) = airefi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    220   cufi(1:klon_omp) = cufi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    221   cvfi(1:klon_omp) = cvfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    222   lonfi(1:klon_omp) = lonfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    223   latfi(1:klon_omp) = latfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
    224   boundslonfi(1:klon_omp,:) = boundslonfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
    225   boundslatfi(1:klon_omp,:) = boundslatfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
     59  ! --> initialize physics distribution, global fields and geometry
     60  ! (i.e. things in phy_common or dynphy_lonlat)
     61  CALL inigeomphy(ii,jj,nlayer, &
     62               nbp, communicator, &
     63               rlatu,rlatv, &
     64               rlonu,rlonv, &
     65               aire,cu,cv)
    22666
    227   ! copy over local grid longitudes and latitudes
    228   CALL init_geometry(klon_omp,lonfi,latfi,boundslonfi,boundslatfi, &
    229                      airefi,cufi,cvfi)
    230 
     67  ! --> now initialize things specific to the phymar physics package
     68 
     69!$OMP PARALLEL
    23170
    23271  ! Initialize tracer names, numbers, etc. for physics

  • LMDZ5/branches/testing/libf/phydev/iophy.F90

    r2408 r2594  
    4444                                jj_nb, jj_begin, jj_end, ii_begin, ii_end, &
    4545                                mpi_size, mpi_rank, klon_mpi, &
    46                                 is_sequential, is_south_pole
     46                                is_sequential, is_south_pole_dyn
    4747  USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, klon_glo
    4848  USE print_control_mod, ONLY: lunout, prt_level
    … …  
    141141      write(lunout,*) "init_iophy_new: mpirank=",mpi_rank," data_ibegin=",data_ibegin," data_iend=",data_iend
    142142      write(lunout,*) "init_iophy_new: mpirank=",mpi_rank," data_ibegin=",data_ibegin," data_iend=",data_iend
    143       write(lunout,*) "init_iophy_new: mpirank=",mpi_rank," is_south_pole=",is_south_pole
     143      write(lunout,*) "init_iophy_new: mpirank=",mpi_rank," is_south_pole=",is_south_pole_dyn
    144144    endif
    145145
    … …  
    148148                            1, nbp_lon, ii_begin, ii_end, jj_begin, jj_end,             &
    149149                            klon_mpi+2*(nbp_lon-1), data_ibegin, data_iend,             &
    150                             io_lat, io_lon,is_south_pole,mpi_rank)
     150                            io_lat, io_lon,is_south_pole_dyn,mpi_rank)
    151151#endif
    152152!$OMP END MASTER

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

    r2435 r2594  
    8888
    8989! Number of diagnostics wavelengths (5 SW + 1 LW @ 10 um)
    90   INTEGER, PARAMETER :: nwave = 5
     90  INTEGER, PARAMETER :: nwave_sw = 5
    9191  INTEGER, PARAMETER :: nwave_lw = 1
     92  INTEGER, PARAMETER :: nwave = nwave_sw + nwave_lw
    9293
    9394! Number of modes spectral bands

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

    r2408 r2594  
    11281128  ENDDO
    11291129
    1130   inu=1
     1130  inu=1         ! visible wavaband
     1131  mrfspecies=2  ! total aerosol AER     
    11311132  DO i=1, KLON
    1132      absvisaer(i)=SUM((1-piz_allaer(i,:,:,inu))*tau_allaer(i,:,:,inu))
     1133     absvisaer(i)=SUM((1-piz_allaer(i,:,mrfspecies,inu))*tau_allaer(i,:,mrfspecies,inu))
    11331134  ENDDO
     1135  STOP
    11341136
    11351137  DEALLOCATE(aerosol_name)

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

    r2408 r2594  
    7474  ! Local
    7575  !
    76   INTEGER, PARAMETER :: las = nwave
     76  INTEGER, PARAMETER :: las = nwave_sw
    7777  LOGICAL :: soluble
    7878 

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

    r2542 r2594  
    44                    wake_pe, wake_fip,  &
    55                    Ale_bl, Ale_bl_trig, Alp_bl, &
    6                     Ale, Alp )
     6                    Ale, Alp, Ale_wake, Alp_wake )
    77
    88! **************************************************************
    … …  
    4646!----------------
    4747  REAL, DIMENSION(klon), INTENT(OUT)                         :: Ale, Alp
     48  REAL, DIMENSION(klon), INTENT(OUT)                         :: Ale_wake, Alp_wake
    4849
    4950  include "thermcell.h"
    … …  
    5556  INTEGER                                                    :: i, k
    5657  REAL, DIMENSION(klon)                                      :: www
    57   REAL, DIMENSION(klon)                                      :: ale_wake, alp_wake
    5858  REAL, SAVE                                                 :: ale_max=1000.
    5959  REAL, SAVE                                                 :: alp_max=2.

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

    r2542 r2594  
    1 SUBROUTINE alpale_th ( dtime, lmax_th, t_seri, &
     1SUBROUTINE alpale_th ( dtime, lmax_th, t_seri, cell_area, &
    22                       cin, s2, n2,  &
    33                       ale_bl_trig, ale_bl_stat, ale_bl,  &
    4                        alp_bl, alp_bl_stat )
     4                       alp_bl, alp_bl_stat, &
     5                       proba_notrig, random_notrig)
    56
    67! **************************************************************
    … …  
    2728!----------------
    2829  REAL, INTENT(IN)                                           :: dtime
     30  REAL, DIMENSION(klon), INTENT(IN)                          :: cell_area
    2931  INTEGER, DIMENSION(klon), INTENT(IN)                       :: lmax_th
    3032  REAL, DIMENSION(klon,klev), INTENT(IN)                     :: t_seri
    … …  
    3739  REAL, DIMENSION(klon), INTENT(INOUT)                       :: alp_bl_stat
    3840
     41  REAL, DIMENSION(klon), INTENT(OUT)                         :: proba_notrig, random_notrig
     42
    3943  include "thermcell.h"
    4044
    … …  
    4448  LOGICAL, SAVE                                              :: first = .TRUE.
    4549  REAL, SAVE                                                 :: random_notrig_max=1.
    46   REAL, DIMENSION(klon)                                      :: proba_notrig, tau_trig, random_notrig
     50  REAL, SAVE                                                 :: cv_feed_area
     51  REAL                                                       :: birth_number
     52  REAL, DIMENSION(klon)                                      :: ale_bl_ref
     53  REAL, DIMENSION(klon)                                      :: tau_trig
     54  REAL, DIMENSION(klon)                                      :: birth_rate
    4755!
    4856    !$OMP THREADPRIVATE(random_notrig_max)
     57    !$OMP THREADPRIVATE(cv_feed_area)
    4958    !$OMP THREADPRIVATE(first)
    5059!
     60 REAL umexp  ! expression of (1.-exp(-x))/x valid for all x, especially when x->0
     61 REAL x
     62 umexp(x) = max(sign(1.,x-1.e-3),0.)*(1.-exp(-x))/max(x,1.e-3) + &
     63            (1.-max(sign(1.,x-1.e-3),0.))*(1.-0.5*x*(1.-x/3.*(1.-0.25*x)))
     64!
     65!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     66!  JYG, 20160513 : Introduction of the Effective Lifting Power (ELP), which
     67! takes into account the area (cv_feed_area) covered by thermals contributing
     68! to each cumulonimbus.
     69!   The use of ELP prevents singularities when the trigger probability tends to
     70! zero. It is activated by iflag_clos_bl = 3.
     71!   The ELP values are stored in the ALP_bl variable.
     72!   
     73!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     74!
     75!---------------------------------------
     76  IF (iflag_clos_bl .LT. 3) THEN
     77!---------------------------------------
     78!
     79!      Original code (Nicolas Rochetin)
     80!     --------------------------------
     81
    5182    IF (first) THEN
     83       random_notrig_max=1.
    5284       CALL getin_p('random_notrig_max',random_notrig_max)
    5385       first=.FALSE.
    … …  
    155187          endif !(iflag_clos_bl)
    156188
     189!
     190!---------------------------------------
     191  ELSEIF (iflag_clos_bl .EQ. 3) THEN  ! (iflag_clos_bl .LT. 3)
     192!---------------------------------------
     193!
     194!      New code with Effective Lifting Power
     195!     -------------------------------------
     196    IF (first) THEN
     197       cv_feed_area = 1.e10   ! m2
     198       CALL getin_p('cv_feed_area', cv_feed_area)
     199       first=.FALSE.
     200    ENDIF
     201
     202          !-----------Stochastic triggering-----------
     203     if (iflag_trig_bl.ge.1) then
     204        !
     205        IF (prt_level .GE. 10) THEN
     206           print *,'cin, ale_bl_stat, alp_bl_stat ', &
     207                cin, ale_bl_stat, alp_bl_stat
     208        ENDIF
     209
     210        ! Use ale_bl_stat (Rochetin's code) or ale_bl (old code) according to
     211        ! iflag_trig_bl value.
     212        IF (iflag_trig_bl.eq.1) then         ! use ale_bl_stat (Rochetin computation)
     213         do i=1,klon
     214              ale_bl_ref(i)=ale_bl_stat(i)
     215         enddo
     216        ELSE IF (iflag_trig_bl.ge.2) then    ! use ale_bl (old computation)
     217         do i=1,klon
     218              ale_bl_ref(i)=Ale_bl(i)
     219         enddo
     220        ENDIF ! (iflag_trig_bl.eq.1)
     221
     222
     223        !----Initializations and random number generation
     224        do i=1,klon
     225           proba_notrig(i)=1.
     226           random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
     227           if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
     228              tau_trig(i)=tau_trig_shallow
     229           else
     230              tau_trig(i)=tau_trig_deep
     231           endif
     232        enddo
     233        !
     234        IF (prt_level .GE. 10) THEN
     235           print *,'random_notrig, tau_trig ', &
     236                random_notrig, tau_trig
     237           print *,'s_trig,s2,n2 ', &
     238                s_trig,s2,n2
     239        ENDIF
     240
     241        !----alp_bl computation
     242        do i=1,klon
     243           if ( (ale_bl_ref(i) .gt. abs(cin(i))+1.e-10) )  then
     244              birth_number = n2(i)*exp(-s_trig/s2(i))
     245              birth_rate(i) = birth_number/(tau_trig(i)*cell_area(i))
     246              proba_notrig(i)=exp(-birth_number*dtime/tau_trig(i))
     247              Alp_bl(i) = Alp_bl(i)* &
     248                          umexp(-birth_number*cv_feed_area/cell_area(i))/ &
     249                          umexp(-birth_number*dtime/tau_trig(i))*  &
     250                          tau_trig(i)*cv_feed_area/(dtime*cell_area(i))
     251          else
     252              proba_notrig(i)=1.
     253              random_notrig(i)=0.
     254              alp_bl(i)=0.
     255           endif
     256        enddo
     257
     258        !----ale_bl_trig computation
     259         do i=1,klon
     260           if (random_notrig(i) .ge. proba_notrig(i)) then
     261              ale_bl_trig(i)=ale_bl_ref(i)
     262           else
     263              ale_bl_trig(i)=0.
     264           endif
     265         enddo
     266
     267        !
     268        IF (prt_level .GE. 10) THEN
     269           print *,'proba_notrig, ale_bl_trig ', &
     270                proba_notrig, ale_bl_trig
     271        ENDIF
     272
     273     endif !(iflag_trig_bl .ge. 1)
     274
     275!---------------------------------------
     276  ENDIF ! (iflag_clos_bl .LT. 3)
     277!---------------------------------------
     278
    157279          IF (prt_level .GE. 10) THEN
    158280             print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
    … …  
    160282
    161283          !cc fin nrlmd le 10/04/2012
    162 
    163           ! ------------------------------------------------------------------
    164           ! Transport de la TKE par les panaches thermiques.
    165           ! FH : 2010/02/01
    166           !     if (iflag_pbl.eq.10) then
    167           !     call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
    168           !    s           rg,paprs,pbl_tke)
    169           !     endif
    170           ! -------------------------------------------------------------------
     284!
    171285          !IM/FH: 2011/02/23
    172286          ! Couplage Thermiques/Emanuel seulement si T<0
    … …  
    180294                endif
    181295             enddo
     296    print *,'In order to run with iflag_coupl=2, you have to comment out the following stop'
     297             STOP
    182298          endif
    183 !
    184299   RETURN
    185300   END

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

    r2542 r2594  
    4040!IM, MAFo fmagic, pmagic : parametres - additionnel et multiplicatif -
    4141!                          pour regler l albedo sur ocean
     42       REAL pbl_lmixmin_alpha
    4243       REAL fmagic, pmagic
    4344! Hauteur (imposee) du contenu en eau du sol
    … …  
    7778       REAL ecrit_LES
    7879       REAL freq_ISCCP, ecrit_ISCCP
    79        REAL freq_COSP
     80       REAL freq_COSP, freq_AIRS
    8081       LOGICAL :: ok_cosp,ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP
     82       LOGICAL :: ok_airs
    8183       INTEGER :: ip_ebil_phy, iflag_rrtm, iflag_ice_thermo, NSW, iflag_albedo
    8284       LOGICAL :: ok_chlorophyll
    … …  
    9496     &     , CH4_ppb, N2O_ppb, CFC11_ppt, CFC12_ppt                     &
    9597     &     , CH4_ppb_per, N2O_ppb_per, CFC11_ppt_per, CFC12_ppt_per     &
    96      &     , cdmmax, cdhmax, ksta, ksta_ter, f_ri_cd_min                &
     98     &     , cdmmax,cdhmax,ksta,ksta_ter,f_ri_cd_min,pbl_lmixmin_alpha  &
    9799     &     , fmagic, pmagic                                             &
    98100     &     , f_cdrag_ter,f_cdrag_oce,f_rugoro,z0min                     &
    … …  
    101103     &     , pasphys            , freq_outNMC, freq_calNMC              &
    102104     &     , lonmin_ins, lonmax_ins, latmin_ins, latmax_ins             &
    103      &     , freq_ISCCP, ecrit_ISCCP, freq_COSP                         &
     105     &     , freq_ISCCP, ecrit_ISCCP, freq_COSP, freq_AIRS              &
    104106     &     , cvl_corr                                                   &
    105107     &     , qsol0,albsno0,evap0                                        &
    … …  
    120122     &     , lev_histins, lev_histLES, lev_histdayNMC, levout_histNMC   &
    121123     &     , ok_histNMC                                                 &
    122      &     , type_run, ok_regdyn, ok_cosp                               &
     124     &     , type_run, ok_regdyn, ok_cosp, ok_airs                      &
    123125     &     , ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP                     &
    124126     &     , ip_ebil_phy                                                &

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

    r2544 r2594  
    55
    66
    7       USE IOIPSL, ONLY : getin
    87      IMPLICIT NONE
    98
    … …  
    2019#include "FCTTRE.h"
    2120#include "thermcell.h"
     21#include "nuage.h"
    2222
    2323      INTEGER itap,ind1,ind2
    … …  
    5454      REAL sth,senv,sigma1s,sigma2s,xth,xenv
    5555      REAL Tbef,zdelta,qsatbef,zcor
    56       REAL alpha,qlbef 
     56      REAL qlbef 
    5757      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
    5858     
    … …  
    6262      REAL erf
    6363
    64       REAL, SAVE :: iflag_cloudth_vert, iflag_cloudth_vert_omp=0
    65 
    66 
    67       LOGICAL, SAVE :: first=.true.
    68 
    69 
    7064
    7165
    7266
    7367!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    74 ! Astuce pour gérer deux versions de cloudth en attendant
    75 ! de converger sur une version nouvelle
     68! Gestion de deux versions de cloudth
    7669!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    77       IF (first) THEN
    78      !$OMP MASTER
    79      CALL getin('iflag_cloudth_vert',iflag_cloudth_vert_omp)
    80      !$OMP END MASTER
    81      !$OMP BARRIER
    82      iflag_cloudth_vert=iflag_cloudth_vert_omp
    83       first=.false.
    84       ENDIF
    85        IF (iflag_cloudth_vert==1) THEN
    86        CALL cloudth_vert(ngrid,klev,ind2,  &
     70
     71      IF (iflag_cloudth_vert.GE.1) THEN
     72      CALL cloudth_vert(ngrid,klev,ind2,  &
    8773     &           ztv,po,zqta,fraca, &
    8874     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
    8975     &           ratqs,zqs,t)
    90        RETURN
    91        ENDIF
     76      RETURN
     77      ENDIF
    9278!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    93 
    9479
    9580
    … …  
    185170      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    186171      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
    187 !      ctot(ind1,ind2)=alpha*cth(ind1,ind2)+(1.-1.*alpha)*cenv(ind1,ind2)
    188 
    189 
    190172
    191173      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
    192174      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
    193175      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    194 !      qltot(ind1,ind2)=alpha*qlth(ind1,ind2)+(1.-1.*alpha)*qlenv(ind1,ind2)
    195      
    196 
    197 !      print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'
    198176
    199177!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    269247     &           ratqs,zqs,t)
    270248
    271 
    272       IMPLICIT NONE
    273 
    274 
    275249!===========================================================================
    276250! Auteur : Arnaud Octavio Jam (LMD/CNRS)
    … …  
    280254
    281255
     256      USE ioipsl_getin_p_mod, ONLY : getin_p
     257
     258      IMPLICIT NONE
     259
    282260#include "YOMCST.h"
    283261#include "YOETHF.h"
    284262#include "FCTTRE.h"
    285263#include "thermcell.h"
    286 
     264#include "nuage.h"
     265     
    287266      INTEGER itap,ind1,ind2
    288267      INTEGER ngrid,klev,klon,l,ig
    … …  
    320299      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
    321300      REAL Tbef,zdelta,qsatbef,zcor
    322       REAL alpha,qlbef 
     301      REAL qlbef 
    323302      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
     303      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
     304      ! (J Jouhaud, JL Dufresne, JB Madeleine)
     305      REAL,SAVE :: vert_alpha
     306      LOGICAL, SAVE :: firstcall = .TRUE.
    324307     
    325308      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
    … …  
    327310      REAL zqs(ngrid), qcloud(ngrid)
    328311      REAL erf
    329 
    330 
    331 
    332 
    333312
    334313!------------------------------------------------------------------------------
    … …  
    355334      sqrtpi=sqrt(pi)
    356335
    357 
     336      IF (firstcall) THEN
     337        vert_alpha=0.5
     338        CALL getin_p('cloudth_vert_alpha',vert_alpha)
     339        WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
     340        firstcall=.FALSE.
     341      ENDIF
    358342
    359343!-------------------------------------------------------------------------------
    … …  
    425409      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    426410      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
    427 !      ctot(ind1,ind2)=alpha*cth(ind1,ind2)+(1.-1.*alpha)*cenv(ind1,ind2)
    428 
    429 
    430411
    431412      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
    432413      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
    433414      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    434 !      qltot(ind1,ind2)=alpha*qlth(ind1,ind2)+(1.-1.*alpha)*qlenv(ind1,ind2)
    435415     
    436 
    437 !      print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'
    438 
    439 
     416       IF (iflag_cloudth_vert == 1) THEN
    440417!-------------------------------------------------------------------------------
    441418!  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
    … …  
    479456      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    480457
     458      ELSE IF (iflag_cloudth_vert == 2) THEN
     459
     460!-------------------------------------------------------------------------------
     461!  Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
     462!-------------------------------------------------------------------------------
     463!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
     464!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
     465!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
     466!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
     467      deltasenv=aenv*vert_alpha*sigma1s
     468      deltasth=ath*vert_alpha*sigma2s
     469     
     470      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
     471      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
     472      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
     473      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
     474!     coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
     475!     coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
     476     
     477      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
     478      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
     479      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
     480
     481      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2)
     482      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
     483      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
     484      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2))
     485
     486!      IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
     487!      IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
     488!      IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
     489
     490      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
     491!      qlenv(ind1,ind2)=IntJ
     492!      print*, qlenv(ind1,ind2),'VERIF EAU'
     493
     494      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2)
     495      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
     496      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
     497      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2))
     498     
     499      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
     500!      qlth(ind1,ind2)=IntJ
     501!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
     502      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
     503     
     504
     505
     506
     507      ENDIF ! of if (iflag_cloudth_vert==1 or 2)
     508
    481509!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     510
    482511      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
    483512      ctot(ind1,ind2)=0.

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

    r2408 r2594  
    164164               iflag_pbl)
    165165       ELSE IF (iflag_pbl<20) THEN
    166           CALL yamada4(knon,dtime,RG,RD,ypaprs,yt, &
     166          CALL yamada4(ni,nsrf,knon,dtime,RG,RD,ypaprs,yt, &
    167167               yzlev,yzlay,yu,yv,yteta, &
    168168               ycdragm,yq2,ykmm,ykmn,ykmq,yustar, &

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

    r2544 r2594  
    9898    REAL,SAVE           :: bl95_b0_omp, bl95_b1_omp
    9999    REAL,SAVE           :: freq_ISCCP_omp, ecrit_ISCCP_omp
    100     REAL,SAVE           :: freq_COSP_omp
     100    REAL,SAVE           :: freq_COSP_omp, freq_AIRS_omp
    101101    real,SAVE           :: fact_cldcon_omp, facttemps_omp,ratqsbas_omp
    102102    real,SAVE           :: tau_cld_cv_omp, coefw_cld_cv_omp
    … …  
    165165    INTEGER,SAVE :: iflag_ice_thermo_omp
    166166    INTEGER,SAVE :: iflag_t_glace_omp
     167    INTEGER,SAVE :: iflag_cloudth_vert_omp
    167168    REAL,SAVE :: rad_froid_omp, rad_chau1_omp, rad_chau2_omp
    168169    REAL,SAVE :: t_glace_min_omp, t_glace_max_omp
    … …  
    178179    REAL,SAVE :: cdmmax_omp,cdhmax_omp,ksta_omp,ksta_ter_omp,f_ri_cd_min_omp
    179180    LOGICAL,SAVE :: ok_kzmin_omp
     181    REAL, SAVE   :: pbl_lmixmin_alpha_omp
    180182    REAL, SAVE ::  fmagic_omp, pmagic_omp
    181183    INTEGER,SAVE :: iflag_pbl_omp,lev_histhf_omp,lev_histday_omp,lev_histmth_omp
    … …  
    188190    REAL, SAVE :: freq_outNMC_omp(3), freq_calNMC_omp(3)
    189191    CHARACTER*4, SAVE :: type_run_omp
    190     LOGICAL,SAVE :: ok_cosp_omp
     192    LOGICAL,SAVE :: ok_cosp_omp, ok_airs_omp
    191193    LOGICAL,SAVE :: ok_mensuelCOSP_omp,ok_journeCOSP_omp,ok_hfCOSP_omp
    192194    REAL,SAVE :: lonmin_ins_omp, lonmax_ins_omp, latmin_ins_omp, latmax_ins_omp
    … …  
    460462    call getin('freq_COSP', freq_COSP_omp)
    461463
     464    !Config Key  = freq_AIRS
     465    !Config Desc = Frequence d'appel du simulateur AIRS en secondes;
     466    !              par defaut 10800, i.e. 3 heures
     467    !Config Def  = 10800.
     468    !Config Help = Used in ini_histdayAIRS.h
     469    !
     470    freq_AIRS_omp = 10800.
     471    call getin('freq_AIRS', freq_AIRS_omp)
     472
    462473    !
    463474    !Config Key  = ip_ebil_phy
    … …  
    11611172
    11621173    !
     1174    !Config Key  = iflag_cloudth_vert
     1175    !Config Desc = 
     1176    !Config Def  = 0
     1177    !Config Help =
     1178    !
     1179    iflag_cloudth_vert_omp = 0
     1180    call getin('iflag_cloudth_vert',iflag_cloudth_vert_omp)
     1181
     1182    !
    11631183    !Config Key  = iflag_ice_thermo
    11641184    !Config Desc = 
    … …  
    12591279    ok_kzmin_omp = .true.
    12601280    call getin('ok_kzmin',ok_kzmin_omp)
     1281
     1282    pbl_lmixmin_alpha_omp=0.0
     1283    call getin('pbl_lmixmin_alpha',pbl_lmixmin_alpha_omp)
     1284
    12611285
    12621286    !
    … …  
    15941618    ok_cosp_omp = .false.
    15951619    call getin('ok_cosp',ok_cosp_omp)
     1620
     1621    !
     1622    !Config Key  = ok_airs
     1623    !Config Desc =
     1624    !Config Def  = .false.
     1625    !Config Help =
     1626    !
     1627    ok_airs_omp = .false.
     1628    call getin('ok_airs',ok_airs_omp)
    15961629
    15971630    !
    … …  
    20442077    exposant_glace = exposant_glace_omp
    20452078    iflag_t_glace = iflag_t_glace_omp
     2079    iflag_cloudth_vert=iflag_cloudth_vert_omp
    20462080    iflag_ice_thermo = iflag_ice_thermo_omp
    20472081    rei_min = rei_min_omp
    … …  
    20552089    f_ri_cd_min = f_ri_cd_min_omp
    20562090    ok_kzmin = ok_kzmin_omp
     2091    pbl_lmixmin_alpha=pbl_lmixmin_alpha_omp
    20572092    fmagic = fmagic_omp
    20582093    pmagic = pmagic_omp
    … …  
    20942129    ecrit_ISCCP = ecrit_ISCCP_omp
    20952130    freq_COSP = freq_COSP_omp
     2131    freq_AIRS = freq_AIRS_omp
    20962132    ok_ade = ok_ade_omp
    20972133    ok_aie = ok_aie_omp
    … …  
    21442180    type_run = type_run_omp
    21452181    ok_cosp = ok_cosp_omp
     2182    ok_airs = ok_airs_omp
     2183
    21462184    ok_mensuelCOSP = ok_mensuelCOSP_omp
    21472185    ok_journeCOSP = ok_journeCOSP_omp
    … …  
    22912329    write(lunout,*)' Frequence appel simulateur ISCCP, ecrit_ISCCP =', ecrit_ISCCP
    22922330    write(lunout,*)' Frequence appel simulateur COSP, freq_COSP =', freq_COSP
     2331    write(lunout,*)' Frequence appel simulateur AIRS, freq_AIRS =', freq_AIRS
    22932332    write(lunout,*)' Sortie bilan d''energie, ip_ebil_phy =', ip_ebil_phy
    22942333    write(lunout,*)' Excentricite = ',R_ecc
    … …  
    23612400    write(lunout,*)' exposant_glace = ',exposant_glace
    23622401    write(lunout,*)' iflag_t_glace = ',iflag_t_glace
     2402    write(lunout,*)' iflag_cloudth_vert = ',iflag_cloudth_vert
    23632403    write(lunout,*)' iflag_ice_thermo = ',iflag_ice_thermo
    23642404    write(lunout,*)' rei_min = ',rei_min
    … …  
    23712411    write(lunout,*)' f_ri_cd_min = ',f_ri_cd_min
    23722412    write(lunout,*)' ok_kzmin = ',ok_kzmin
     2413    write(lunout,*)' pbl_lmixmin_alpha = ',pbl_lmixmin_alpha
    23732414    write(lunout,*)' fmagic = ',fmagic
    23742415    write(lunout,*)' pmagic = ',pmagic
    … …  
    24042445    write(lunout,*)' type_run = ',type_run
    24052446    write(lunout,*)' ok_cosp = ',ok_cosp
     2447    write(lunout,*)' ok_airs = ',ok_airs
     2448
    24062449    write(lunout,*)' ok_mensuelCOSP = ',ok_mensuelCOSP
    24072450    write(lunout,*)' ok_journeCOSP = ',ok_journeCOSP

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

    r2435 r2594  
    108108                   -31.5,-28.5,-25.5,-22.5,-19.5,-16.5,-13.5,-10.5, -7.5, -4.5, &
    109109                    -1.5,  1.5,  4.5,  7.5, 10.5, 13.5, 16.5, 19.5, 22.5, 25.5/)
    110     real,parameter,dimension(2,LIDAR_NTEMP) :: LIDAR_PHASE_TEMP_BNDS=reshape(source=(/-273.15,-90.,-90.,-87.,-87.,-84.,-84.,-81.,-81.,-78., &
     110    real,parameter,dimension(2,LIDAR_NTEMP) :: LIDAR_PHASE_TEMP_BNDS=reshape(source=&
     111                    (/-273.15,-90.,-90.,-87.,-87.,-84.,-84.,-81.,-81.,-78., &
    111112                   -78.,-75.,-75.,-72.,-72.,-69.,-69.,-66.,-66.,-63., &
    112113                   -63.,-60.,-60.,-57.,-57.,-54.,-54.,-51.,-51.,-48., &

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

    r2471 r2594  
    228228  real,dimension(2,SR_BINS) :: sratio_bounds
    229229  real,dimension(SR_BINS)   ::  sratio_ax
    230   CHARACTER(LEN=20), DIMENSION(3)  :: chfreq = (/ '1day', '1d', '3h' /)           
     230  CHARACTER(LEN=20), DIMENSION(3)  :: chfreq = (/ '1day', '1d  ', '3h  ' /)           
    231231
    232232!!! Variables d'entree
    … …  
    363363      CALL histvert(cosp_nidfiles(iff),"column","column","count",Ncolumns,column_ax(1:Ncolumns),nvertcol(iff))
    364364
    365       CALL histvert(cosp_nidfiles(iff),"temp","temperature","C",LIDAR_NTEMP,LIDAR_PHASE_TEMP,nverttemp(iff))                                       
    366       CALL histvert(cosp_nidfiles(iff),"cth16","altitude","m",MISR_N_CTH,MISR_CTH,nvertmisr(iff))                                                                                                 
     365      CALL histvert(cosp_nidfiles(iff),"temp","temperature","C",LIDAR_NTEMP,LIDAR_PHASE_TEMP,nverttemp(iff))
     366
     367      CALL histvert(cosp_nidfiles(iff),"cth16","altitude","m",MISR_N_CTH,MISR_CTH,nvertmisr(iff))
     368
    367369!      CALL histvert(cosp_nidfiles(iff),"dbze","equivalent_reflectivity_factor","dBZ",DBZE_BINS,dbze_ax,nvertbze(iff))
    368370     

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

    r2488 r2594  
    322322    if(modis%Cloud_Particle_Size_Water_Mean(ip).eq.R_UNDEF)then
    323323       modis%Cloud_Particle_Size_Water_Mean(ip)=Cosp_fill_value
     324    endif
     325    if(modis%Cloud_Particle_Size_Ice_Mean(ip).eq.R_UNDEF)then
     326       modis%Cloud_Particle_Size_Ice_Mean(ip)=Cosp_fill_value
    324327    endif
    325328    if(modis%Cloud_Top_Pressure_Total_Mean(ip).eq.R_UNDEF)then

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

    r2435 r2594  
    7373
    7474!! AI rajouter
    75   #include "cosp_defs.h"
     75#include "cosp_defs.h"
    7676  USE MOD_COSP_CONSTANTS
    7777  USE MOD_COSP_TYPES

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

    r2542 r2594  
    38603860      sigd(il) = sigd(il)/alpha_qpos(il)
    38613861      precip(il) = precip(il)/alpha_qpos(il)
     3862      cbmf(il) = cbmf(il)/alpha_qpos(il)
    38623863    END IF
    38633864  END DO

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

    r2471 r2594  
    32363236     &           thlpcar(nlev_max),tracer(nlev_max,ntrac)
    32373237
     3238        real height1(nlev_max)
     3239
    32383240        integer, parameter :: ilesfile=1
    32393241        integer :: ierr,k,itrac,nt1,nt2
    … …  
    32453247        read (ilesfile,*) kmax
    32463248        do k=1,kmax
    3247           read (ilesfile,*) height(k),thlprof(k),qtprof (k),               &
     3249          read (ilesfile,*) height1(k),thlprof(k),qtprof (k),               &
    32483250     &                      uprof (k),vprof  (k),e12prof(k)
    32493251        enddo
    … …  
    32613263          read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k),         &
    32623264     &                      dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k)
     3265        end do
     3266        do k=1,kmax
     3267          if (height(k) .ne. height1(k)) then
     3268            print *, 'fichiers prof.inp et lscale.inp incompatibles :'
     3269            print *, 'les niveaux different : ',k,height1(k), height(k)
     3270            stop
     3271          endif
    32633272        end do
    32643273        close(ilesfile)

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

    r2408 r2594  
    22! $Id$
    33!
     4!---------------------------------------------------------------------
     5! Forcing_LES case: constant dq_dyn
     6!---------------------------------------------------------------------
     7      if (forcing_LES) then
     8        DO l = 1,llm
     9          d_q_adv(l,1) = dq_dyn(l,1)
     10        ENDDO
     11      endif ! forcing_LES
     12!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    413!---------------------------------------------------------------------
    514! Interpolation forcing in time and onto model levels

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

    r2488 r2594  
    3333   USE mod_1D_cases_read
    3434   USE mod_1D_amma_read
    35    USE print_control_mod, ONLY: prt_level
     35   USE print_control_mod, ONLY: lunout, prt_level
    3636   USE iniphysiq_mod, ONLY: iniphysiq
    3737   USE mod_const_mpi, ONLY: comm_lmdz
    … …  
    5454#include "fcg_gcssold.h"
    5555!!!#include "fbforcing.h"
     56#include "compbl.h"
    5657
    5758!=====================================================================
    … …  
    242243!---------------------------------------------------------------------
    243244      integer :: k,l,i,it=1,mxcalc
     245      integer :: nsrf
    244246      integer jcode
    245247      INTEGER read_climoz
    … …  
    738740        t_ancien(1,:)=temp(:)
    739741        q_ancien(1,:)=q(:,1)
    740         pbl_tke(:,:,:)=1.e-8
    741         wake_delta_pbl_tke(:,:,:)=0.
     742!jyg<
     743!!        pbl_tke(:,:,:)=1.e-8
     744        pbl_tke(:,:,:)=0.
     745        pbl_tke(:,2,:)=1.e-2
     746        PRINT *, ' pbl_tke dans lmdz1d '
     747       DO nsrf = 1,4
     748         PRINT *,'pbl_tke(1,:,',nsrf,') ',pbl_tke(1,:,nsrf)
     749       ENDDO
     750
     751!>jyg
    742752
    743753        rain_fall=0.
    … …  
    762772        wake_deltaq = 0.
    763773        wake_deltat = 0.
    764         wake_delta_pbl_TKE = 0.
     774        wake_delta_pbl_TKE(:,:,:) = 0.
    765775        delta_tsurf = 0.
    766776        wake_fip = 0.
    … …  
    791801! run_off_lic_0,pbl_tke(:,1:klev,nsrf), zmax0,f0,sig1,w01
    792802! wake_deltat,wake_deltaq,wake_s,wake_cstar,wake_fip,wake_delta_pbl_tke(:,1:klev,nsrf)
     803!
     804! NB2: The content of the startphy.nc file depends on some flags defined in
     805! the ".def" files. However, since conf_phys is not called in lmdz1d.F90, these flags have
     806! to be set at some arbitratry convenient values.
    793807!------------------------------------------------------------------------
    794808!Al1 =============== restart option ==========================
    795809        if (.not.restart) then
     810          iflag_pbl = 5
    796811          call phyredem ("startphy.nc")
    797812        else
    … …  
    10751090!---------------------------------------------------------------------
    10761091
    1077       IF (nudge_tsoil) THEN
     1092      IF (nudge_tsoil .AND. .NOT. lastcall) THEN
    10781093       ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:)                     &
    10791094     &  -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)

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

    r2408 r2594  
    33!*******************************************************************************
    44
     5  USE print_control_mod, ONLY: lunout
     6  USE assert_eq_m,       ONLY: assert_eq
    57  PRIVATE
    6   PUBLIC :: grid_noro
     8  PUBLIC :: grid_noro, grid_noro0
    79
    810
    … …  
    4547!    on the edge are contained in input cells.
    4648!===============================================================================
    47   USE assert_eq_m, ONLY: assert_eq
    48   USE print_control_mod, ONLY: lunout
    49   IMPLICIT NONE
    50   REAL, PARAMETER :: epsfra = 1.e-5
     49  IMPLICIT NONE
    5150!-------------------------------------------------------------------------------
    5251! Arguments:
    … …  
    315314!-------------------------------------------------------------------------------
    316315!
     316SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)
     317!
     318!===============================================================================
     319! Purpose: Extracted from grid_noro to provide geopotential height for dynamics
     320!          without any call to physics subroutines.
     321!===============================================================================
     322  IMPLICIT NONE
     323!-------------------------------------------------------------------------------
     324! Arguments:
     325  REAL, INTENT(IN)   :: xd(:), yd(:) !--- INPUT  COORDINATES     (imdp) (jmdp)
     326  REAL, INTENT(IN)   :: zd(:,:)      !--- INPUT  FIELD           (imdp,jmdp)
     327  REAL, INTENT(IN)   :: x(:), y(:)   !--- OUTPUT COORDINATES     (imar+1) (jmar)
     328  REAL, INTENT(OUT)  :: zphi(:,:)    !--- GEOPOTENTIAL           (imar+1,jmar)
     329  REAL, INTENT(INOUT):: mask(:,:)    !--- MASK                   (imar+1,jmar)
     330!-------------------------------------------------------------------------------
     331! Local variables:
     332  CHARACTER(LEN=256) :: modname="grid_noro0"
     333  REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
     334  REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
     335  REAL, ALLOCATABLE :: weight(:,:)                ! dim (imar+1,jmar)
     336  REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:)   ! dim (imar+1,jmar)
     337  REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
     338  REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
     339  REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
     340  LOGICAL :: masque_lu
     341  INTEGER :: i, ii, imdp, imar, iext
     342  INTEGER :: j, jj, jmdp, jmar, nn
     343  REAL    :: xpi, zlenx, weighx, xincr,  zbordnor, zmeanor, zweinor, zbordest
     344  REAL    :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue
     345!-------------------------------------------------------------------------------
     346  imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
     347  jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
     348  imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
     349  jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
     350! IF(imar/=iim)   CALL abort_gcm(TRIM(modname),'imar/=iim'  ,1)
     351! IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
     352  iext=imdp/10
     353  xpi = ACOS(-1.)
     354  rad = 6371229.
     355
     356!--- ARE WE USING A READ MASK ?
     357  masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
     358  WRITE(lunout,*)'Masque lu: ',masque_lu
     359
     360!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
     361  ALLOCATE(xusn(imdp+2*iext))
     362  xusn(1     +iext:imdp  +iext)=xd(:)
     363  xusn(1          :       iext)=xd(1+imdp-iext:imdp)-2.*xpi
     364  xusn(1+imdp+iext:imdp+2*iext)=xd(1          :iext)+2.*xpi
     365
     366  ALLOCATE(yusn(jmdp+2))
     367  yusn(1       )=yd(1)   +(yd(1)   -yd(2))
     368  yusn(2:jmdp+1)=yd(:)
     369  yusn(  jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
     370
     371  ALLOCATE(zusn(imdp+2*iext,jmdp+2))
     372  zusn(1       +iext:imdp  +iext,2:jmdp+1)=zd  (:                   ,     :)
     373  zusn(1            :       iext,2:jmdp+1)=zd  (imdp-iext+1:imdp    ,     :)
     374  zusn(1+imdp  +iext:imdp+2*iext,2:jmdp+1)=zd  (1:iext              ,     :)
     375  zusn(1            :imdp/2+iext,       1)=zusn(1+imdp/2:imdp  +iext,     2)
     376  zusn(1+imdp/2+iext:imdp+2*iext,       1)=zusn(1       :imdp/2+iext,     2)
     377  zusn(1            :imdp/2+iext,  jmdp+2)=zusn(1+imdp/2:imdp  +iext,jmdp+1)
     378  zusn(1+imdp/2+iext:imdp+2*iext,  jmdp+2)=zusn(1       :imdp/2+iext,jmdp+1)
     379
     380!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
     381  ALLOCATE(a(imar+1),b(imar+1))
     382  b(1:imar)=(x(1:imar  )+ x(2:imar+1))/2.0
     383  b(imar+1)= x(  imar+1)+(x(  imar+1)-x(imar))/2.0
     384  a(1)=x(1)-(x(2)-x(1))/2.0
     385  a(2:imar+1)= b(1:imar)
     386
     387  ALLOCATE(c(jmar),d(jmar))
     388  d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
     389  d(  jmar  )= y(  jmar  )+(y(  jmar)-y(jmar-1))/2.0
     390  c(1)=y(1)-(y(2)-y(1))/2.0
     391  c(2:jmar)=d(1:jmar-1)
     392
     393!--- INITIALIZATIONS:
     394  ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
     395  ALLOCATE(zmea  (imar+1,jmar)); zmea  (:,:)= 0.0
     396
     397!--- SUMMATION OVER GRIDPOINT AREA
     398  zleny=xpi/REAL(jmdp)*rad
     399  xincr=xpi/REAL(jmdp)/2.
     400  ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
     401  ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
     402  DO ii = 1, imar+1
     403    DO jj = 1, jmar
     404      DO j = 2,jmdp+1
     405        zlenx  =zleny  *COS(yusn(j))
     406        zbordnor=(xincr+c(jj)-yusn(j))*rad
     407        zbordsud=(xincr-d(jj)+yusn(j))*rad
     408        weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
     409        IF(weighy/=0) THEN
     410          DO i = 2, imdp+2*iext-1
     411            zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))
     412            zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))
     413            weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
     414            IF(weighx/=0)THEN
     415              num_tot(ii,jj)=num_tot(ii,jj)+1.0
     416              IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
     417              weight(ii,jj)=weight(ii,jj)+weighx*weighy
     418              zmea  (ii,jj)=zmea  (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
     419            END IF
     420          END DO
     421        END IF
     422      END DO
     423    END DO
     424  END DO
     425
     426!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
     427  IF(.NOT.masque_lu) THEN
     428    WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
     429  END IF
     430  nn=COUNT(weight(:,1:jmar-1)==0.0)
     431  IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
     432  WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
     433
     434!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)
     435  ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
     436  WHERE(mask>=0.1) mask_tmp = 1.
     437  WHERE(weight(:,:)/=0.0)
     438    zphi(:,:)=mask_tmp(:,:)*zmea(:,:)
     439    zmea(:,:)=mask_tmp(:,:)*zmea(:,:)
     440  END WHERE
     441  WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
     442
     443!--- Values at poles
     444  zphi(imar+1,:)=zphi(1,:)
     445
     446  zweinor=SUM(weight(1:imar,   1),DIM=1)
     447  zweisud=SUM(weight(1:imar,jmar),DIM=1)
     448  zmeanor=SUM(weight(1:imar,   1)*zmea(1:imar,   1),DIM=1)
     449  zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)
     450  zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud
     451
     452END SUBROUTINE grid_noro0
     453!
     454!-------------------------------------------------------------------------------
     455
     456
     457!-------------------------------------------------------------------------------
     458!
    317459SUBROUTINE MVA9(x)
    318460!

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

    r2544 r2594  
    99      REAL tmax_fonte_cv
    1010
    11       INTEGER iflag_t_glace,iflag_cld_cv
     11      INTEGER iflag_t_glace, iflag_cloudth_vert, iflag_cld_cv
    1212
    1313      common /nuagecom/ rad_froid,rad_chau1, rad_chau2,t_glace_max,     &
    1414     &                  t_glace_min,exposant_glace,rei_min,rei_max,     &
    1515     &                  tau_cld_cv,coefw_cld_cv,                        &
    16      &                  iflag_t_glace,iflag_cld_cv,tmax_fonte_cv
     16     &                  iflag_t_glace,iflag_cloudth_vert,iflag_cld_cv,  &
     17     &                  tmax_fonte_cv
    1718!$OMP THREADPRIVATE(/nuagecom/)

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

    r2570 r2594  
    1818       wake_deltat, wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
    1919       wake_s, zgam, zmax0, zmea, zpic, zsig, &
    20        zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl
     20       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, u10m, v10m
    2121  USE geometry_mod, ONLY : longitude_deg, latitude_deg
    2222  USE iostart, ONLY : close_startphy, get_field, get_var, open_startphy
    … …  
    260260  ENDDO
    261261
    262 
    263  
     262  found=phyetat0_srf(1,u10m,"U10M","u a 10m",0.)
     263  found=phyetat0_srf(1,v10m,"V10M","v a 10m",0.)
     264
    264265!===================================================================
    265266  ! Lecture des temperatures du sol profond:

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

    r2570 r2594  
    2121                                wake_pe, wake_fip, fm_therm, entr_therm,     &
    2222                                detr_therm, Ale_bl, Ale_bl_trig, Alp_bl,     &
    23                                 du_gwd_rando, du_gwd_front
     23                                du_gwd_rando, du_gwd_front, u10m, v10m
    2424  USE geometry_mod, ONLY : longitude_deg, latitude_deg
    2525  USE iostart, ONLY: open_restartphy, close_restartphy, put_field, put_var
    … …  
    153153
    154154 
     155  CALL put_field_srf1("U10M", "u a 10m", u10m)
     156
     157  CALL put_field_srf1("V10M", "v a 10m", v10m)
     158
     159
    155160! ================== Tsoil =========================================
    156161  CALL put_field_srf2("Tsoil","Temperature",tsoil(:,:,:))

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

    r2542 r2594  
    1616      REAL, SAVE, ALLOCATABLE :: u_seri(:,:), v_seri(:,:)
    1717      !$OMP THREADPRIVATE(u_seri, v_seri)
     18      REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:), l_mix(:,:,:)
     19      !$OMP THREADPRIVATE(l_mixmin, l_mix)
    1820
    1921      REAL, SAVE, ALLOCATABLE :: tr_seri(:,:,:)
    … …  
    406408      allocate(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev))
    407409      allocate(u_seri(klon,klev),v_seri(klon,klev))
     410      allocate(l_mixmin(klon,klev,nbsrf), l_mix(klon,klev,nbsrf))
     411      l_mix(:,:,:)=0. ; l_mixmin(:,:,:)=0. ! doit etre initialse car pas toujours remplis
    408412
    409413      allocate(tr_seri(klon,klev,nbtr))
    … …  
    625629      deallocate(t_seri,q_seri,ql_seri,qs_seri)
    626630      deallocate(u_seri,v_seri)
     631      deallocate(l_mixmin,l_mix)
    627632
    628633      deallocate(tr_seri)

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

    r2542 r2594  
    4141
    4242!!! 2D
     43
     44! Marine
     45
     46  TYPE(ctrl_out), SAVE :: o_alt_tropo = ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     47  'alt_tropo','Tropopause pressure','hPa',&
     48   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     49    "inst(X)", "inst(X)","inst(X)" /))
     50
     51  TYPE(ctrl_out), SAVE :: o_map_prop_hc = ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     52  'map_prop_hc','Proportion of high clouds',' ',&
     53   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     54    "inst(X)", "inst(X)","inst(X)" /))
     55
     56  TYPE(ctrl_out), SAVE :: o_map_prop_hist = &
     57  ctrl_out((/1,1,1,1,1,1,10,10,10/),&
     58  'map_prop_hist','Proportion of high ice semi-transp clouds',' ',&
     59   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     60    "inst(X)", "inst(X)","inst(X)" /))
     61
     62  TYPE(ctrl_out), SAVE :: o_map_emis_hc = &
     63  ctrl_out((/1,1,1,1,1,1,10,10,10/),&
     64  'map_emis_hc','Emissivity of high clouds',' ',&
     65   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     66    "inst(X)", "inst(X)","inst(X)" /))
     67
     68  TYPE(ctrl_out), SAVE :: o_map_iwp_hc = &
     69  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     70  'map_iwp_hc','Ice water path of high clouds','g/m2',&
     71   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     72    "inst(X)", "inst(X)","inst(X)" /))
     73
     74  TYPE(ctrl_out), SAVE :: o_map_deltaz_hc = &
     75  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     76  'map_deltaz_hc','geom thickness of high clouds','m',&
     77   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     78    "inst(X)", "inst(X)","inst(X)" /))
     79
     80  TYPE(ctrl_out), SAVE :: o_map_pcld_hc = &
     81  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     82  'map_pcld_hc','cloud pressure of high clouds','hPa',&
     83   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     84    "inst(X)", "inst(X)","inst(X)" /))
     85
     86   TYPE(ctrl_out), SAVE :: o_map_tcld_hc = &
     87  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     88  'map_tcld_hc','cloud temperature of high clouds','K',&
     89   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     90    "inst(X)", "inst(X)","inst(X)" /))
     91
     92
     93  TYPE(ctrl_out), SAVE :: o_map_emis_hist = &
     94  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     95  'map_emis_hist','Emissivity of high ice st clouds',' ',&
     96   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     97    "inst(X)", "inst(X)","inst(X)" /))
     98
     99  TYPE(ctrl_out), SAVE :: o_map_iwp_hist = &
     100  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     101  'map_iwp_hist','Ice water path of high ice st clouds','g/m2',&
     102   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     103    "inst(X)", "inst(X)","inst(X)" /))
     104
     105  TYPE(ctrl_out), SAVE :: o_map_deltaz_hist = &
     106  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     107  'map_deltaz_hist','geom thickness of high ice st clouds','m',&
     108   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     109    "inst(X)", "inst(X)","inst(X)" /))
     110
     111  TYPE(ctrl_out), SAVE :: o_map_rad_hist = &
     112  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     113  'map_rad_hist','ice crystals radius in high ice st clouds','µm',&
     114   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     115    "inst(X)", "inst(X)","inst(X)" /))
     116
     117
     118 TYPE(ctrl_out), SAVE :: o_map_emis_Cb = &
     119  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     120  'map_emis_Cb','Emissivity of high Cb clouds',' ',&
     121   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     122    "inst(X)", "inst(X)","inst(X)" /))
     123
     124 TYPE(ctrl_out), SAVE :: o_map_pcld_Cb = &
     125  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     126  'map_pcld_Cb','cloud pressure of high Cb clouds','hPa',&
     127   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     128    "inst(X)", "inst(X)","inst(X)" /))
     129
     130 TYPE(ctrl_out), SAVE :: o_map_tcld_Cb = &
     131  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     132  'map_tcld_Cb','cloud temperature of high Cb clouds','K',&
     133   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     134    "inst(X)", "inst(X)","inst(X)" /))
     135
     136
     137 TYPE(ctrl_out), SAVE :: o_map_emis_Anv = &
     138  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     139  'map_emis_Anv','Emissivity of high Anv clouds',' ',&
     140   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     141    "inst(X)", "inst(X)","inst(X)" /))
     142
     143 TYPE(ctrl_out), SAVE :: o_map_pcld_Anv = &
     144  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     145  'map_pcld_Anv','cloud pressure of high Anv clouds','hPa',&
     146   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     147    "inst(X)", "inst(X)","inst(X)" /))
     148
     149  TYPE(ctrl_out), SAVE :: o_map_tcld_Anv = &
     150  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     151  'map_tcld_Anv','cloud temperature of high Anv clouds','K',&
     152   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     153    "inst(X)", "inst(X)","inst(X)" /))
     154
     155  TYPE(ctrl_out), SAVE :: o_map_emis_ThCi = &
     156  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     157  'map_emis_ThCi','Emissivity of high ThCi clouds',' ',&
     158   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     159    "inst(X)", "inst(X)","inst(X)" /))
     160
     161  TYPE(ctrl_out), SAVE :: o_map_pcld_ThCi = &
     162  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     163  'map_pcld_ThCi','cloud pressure of high ThCi clouds','hPa',&
     164   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     165    "inst(X)", "inst(X)","inst(X)" /))
     166
     167  TYPE(ctrl_out), SAVE :: o_map_tcld_ThCi = &
     168  ctrl_out((/10,10,1,10,10,10,10,10,10/),&
     169  'map_tcld_ThCi','cloud temperature of high ThCi clouds','K',&
     170   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     171    "inst(X)", "inst(X)","inst(X)" /))
     172
     173   TYPE(ctrl_out), SAVE :: o_map_ntot = &
     174  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     175  'map_ntot','total AIRS cloud fraction',' ',&
     176   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     177    "inst(X)", "inst(X)","inst(X)" /))
     178
     179  TYPE(ctrl_out), SAVE :: o_map_hc = &
     180  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     181  'map_hc','high clouds AIRS cloud fraction',' ',&
     182   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     183    "inst(X)", "inst(X)","inst(X)" /))
     184
     185  TYPE(ctrl_out), SAVE :: o_map_hist = &
     186  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     187  'map_hist','high clouds ice st AIRS cloud fraction',' ',&
     188   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     189    "inst(X)", "inst(X)","inst(X)" /))
     190
     191  TYPE(ctrl_out), SAVE :: o_map_Cb = &
     192  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     193  'map_Cb','high clouds Cb AIRS cloud fraction',' ',&
     194   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     195    "inst(X)", "inst(X)","inst(X)" /))
     196
     197 TYPE(ctrl_out), SAVE :: o_map_ThCi = &
     198  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     199  'map_ThCi','high clouds ThCi AIRS cloud fraction',' ',&
     200   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     201    "inst(X)", "inst(X)","inst(X)" /))
     202
     203 TYPE(ctrl_out), SAVE :: o_map_Anv = &
     204  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
     205  'map_Anv','high clouds Anv AIRS cloud fraction',' ',&
     206   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
     207    "inst(X)", "inst(X)","inst(X)" /))
     208
     209
     210! Fin Marine
     211
    43212  TYPE(ctrl_out), SAVE :: o_flat = ctrl_out((/ 5, 1, 10, 10, 5, 10, 11, 11, 11 /), &
    44213    'flat', 'Latent heat flux', 'W/m2', (/ ('', i=1, 9) /))
    … …  
    745914  TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_tke_srf      = (/             &
    746915      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'tke_ter',       &
    747       "Max Turb. Kinetic Energy "//clnsurf(1),"-", (/ ('', i=1, 9) /)), &
     916      "Max Turb. Kinetic Energy "//clnsurf(1),"m2/s2", (/ ('', i=1, 9) /)), &
    748917      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'tke_lic',       &
    749       "Max Turb. Kinetic Energy "//clnsurf(2),"-", (/ ('', i=1, 9) /)), &
     918      "Max Turb. Kinetic Energy "//clnsurf(2),"m2/s2", (/ ('', i=1, 9) /)), &
    750919      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'tke_oce',       &
    751       "Max Turb. Kinetic Energy "//clnsurf(3),"-", (/ ('', i=1, 9) /)), &
     920      "Max Turb. Kinetic Energy "//clnsurf(3),"m2/s2", (/ ('', i=1, 9) /)), &
    752921      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11 /),'tke_sic',       &
    753       "Max Turb. Kinetic Energy "//clnsurf(4),"-", (/ ('', i=1, 9) /)) /)
     922      "Max Turb. Kinetic Energy "//clnsurf(4),"m2/s2", (/ ('', i=1, 9) /)) /)
     923
     924  TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_l_mixmin      = (/             &
     925      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'l_mixmin_ter',       &
     926      "PBL mixing length "//clnsurf(1),"m", (/ ('', i=1, 9) /)), &
     927      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'l_mixmin_lic',       &
     928      "PBL mixing length "//clnsurf(2),"m", (/ ('', i=1, 9) /)), &
     929      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'l_mixmin_oce',       &
     930      "PBL mixing length "//clnsurf(3),"m", (/ ('', i=1, 9) /)), &
     931      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'l_mixmin_sic',       &
     932      "PBL mixing length "//clnsurf(4),"m", (/ ('', i=1, 9) /)) /)
     933
     934  TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_l_mix      = (/             &
     935      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'l_mix_ter',       &
     936      "min PBL mixing length "//clnsurf(1),"m", (/ ('', i=1, 9) /)), &
     937      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'l_mix_lic',       &
     938      "min PBL mixing length "//clnsurf(2),"m", (/ ('', i=1, 9) /)), &
     939      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'l_mix_oce',       &
     940      "min PBL mixing length "//clnsurf(3),"m", (/ ('', i=1, 9) /)), &
     941      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'l_mix_sic',       &
     942      "min PBL mixing length "//clnsurf(4),"m", (/ ('', i=1, 9) /)) /)
    754943
    755944  TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_tke_max_srf  = (/                          &
    … …  
    9681157    'ovapinit', 'Specific humidity (begin of timestep)', 'kg/kg', (/ ('', i=1, 9) /))
    9691158  TYPE(ctrl_out), SAVE :: o_oliq = ctrl_out((/ 2, 3, 4, 10, 10, 10, 11, 11, 11 /), &
    970     'oliq', 'Condensed water', 'kg/kg', (/ ('', i=1, 9) /))
     1159    'oliq', 'Liquid water', 'kg/kg', (/ ('', i=1, 9) /))
     1160  TYPE(ctrl_out), SAVE :: o_ocond = ctrl_out((/ 2, 3, 4, 10, 10, 10, 11, 11, 11 /), &
     1161    'ocond', 'Condensed water', 'kg/kg', (/ ('', i=1, 9) /))
    9711162  TYPE(ctrl_out), SAVE :: o_wvapp = ctrl_out((/ 2, 10, 10, 10, 10, 10, 11, 11, 11 /), &
    9721163    'wvapp', '', '', (/ ('', i=1, 9) /))

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

    r2542 r2594  
    5656    include "clesphys.h"
    5757    include "thermcell.h"
     58    include "YOMCST.h"
    5859
    5960    ! ug Nouveaux arguments n\'ecessaires au histwrite_mod:
    … …  
    263264     ELSE
    264265         CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
    265          CALL ymds2ju(annee_ref, 1, day_ini, start_time, zjulian_start)
     266         CALL ymds2ju(annee_ref, 1, day_ini, start_time*rday, zjulian_start)
    266267     END IF
    267268

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

    r2542 r2594  
    1616  INTEGER, SAVE, ALLOCATABLE ::  itau_con(:)       ! Nombre de pas ou rflag <= 1
    1717  !$OMP THREADPRIVATE(itau_con)
    18   REAL, ALLOCATABLE :: bils_ec(:) ! Contribution of energy conservation
    19   REAL, ALLOCATABLE :: bils_ech(:) ! Contribution of energy conservation
    20   REAL, ALLOCATABLE :: bils_tke(:) ! Contribution of energy conservation
    21   REAL, ALLOCATABLE :: bils_diss(:) ! Contribution of energy conservation
    22   REAL, ALLOCATABLE :: bils_kinetic(:) ! bilan de chaleur au sol, kinetic
    23   REAL, ALLOCATABLE :: bils_enthalp(:) ! bilan de chaleur au sol
    24   REAL, ALLOCATABLE :: bils_latent(:) ! bilan de chaleur au sol
     18  REAL, SAVE, ALLOCATABLE :: bils_ec(:) ! Contribution of energy conservation
     19  REAL, SAVE, ALLOCATABLE :: bils_ech(:) ! Contribution of energy conservation
     20  REAL, SAVE, ALLOCATABLE :: bils_tke(:) ! Contribution of energy conservation
     21  REAL, SAVE, ALLOCATABLE :: bils_diss(:) ! Contribution of energy conservation
     22  REAL, SAVE, ALLOCATABLE :: bils_kinetic(:) ! bilan de chaleur au sol, kinetic
     23  REAL, SAVE, ALLOCATABLE :: bils_enthalp(:) ! bilan de chaleur au sol
     24  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)
    2626
     27! Marine
     28! Variables de sortie du simulateur AIRS
     29
     30  REAL, SAVE, ALLOCATABLE :: map_prop_hc(:),map_prop_hist(:),alt_tropo(:)
     31  !$OMP THREADPRIVATE(map_prop_hc,map_prop_hist,alt_tropo)
     32  REAL, SAVE, ALLOCATABLE :: map_emis_hc(:),map_iwp_hc(:),map_deltaz_hc(:), &
     33                       map_pcld_hc(:),map_tcld_hc(:)
     34  !$OMP THREADPRIVATE(map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc)
     35  REAL, SAVE, ALLOCATABLE :: map_emis_hist(:),map_iwp_hist(:),map_deltaz_hist(:),map_rad_hist(:)         
     36  !$OMP THREADPRIVATE(map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist)
     37  REAL, SAVE, ALLOCATABLE :: map_ntot(:),map_hc(:),map_hist(:)
     38  REAL, SAVE, ALLOCATABLE :: map_Cb(:),map_ThCi(:),map_Anv(:)
     39  !$OMP THREADPRIVATE(map_ntot,map_hc,map_hist,map_Cb,map_ThCi,map_Anv)
     40  REAL, SAVE, ALLOCATABLE :: map_emis_Cb(:),map_pcld_Cb(:),map_tcld_Cb(:)
     41  REAL, SAVE, ALLOCATABLE :: map_emis_ThCi(:),map_pcld_ThCi(:),map_tcld_ThCi(:)
     42  !$OMP THREADPRIVATE(map_emis_Cb,map_pcld_Cb,map_tcld_Cb,map_emis_ThCi)
     43  REAL, SAVE, ALLOCATABLE :: map_emis_Anv(:),map_pcld_Anv(:),map_tcld_Anv(:)
     44  !$OMP THREADPRIVATE(map_pcld_ThCi,map_tcld_ThCi,map_emis_Anv,map_pcld_Anv,map_tcld_Anv)             
     45   
    2746
    2847  ! ug Plein de variables venues de phys_output_mod
    … …  
    103122    allocate (bils_ec(klon),bils_ech(klon),bils_tke(klon),bils_diss(klon),bils_kinetic(klon),bils_enthalp(klon),bils_latent(klon))
    104123
     124! Marine
     125! Variables de sortie simulateur AIRS
     126
     127!     if (ok_airs) then
     128      allocate (map_prop_hc(klon),map_prop_hist(klon))
     129      allocate (alt_tropo(klon))
     130      allocate (map_emis_hc(klon),map_iwp_hc(klon),map_deltaz_hc(klon))
     131      allocate (map_pcld_hc(klon),map_tcld_hc(klon))
     132      allocate (map_emis_hist(klon),map_iwp_hist(klon),map_deltaz_hist(klon))
     133      allocate (map_rad_hist(klon))
     134      allocate (map_ntot(klon),map_hc(klon),map_hist(klon))
     135      allocate (map_Cb(klon),map_ThCi(klon),map_Anv(klon))
     136      allocate (map_emis_Cb(klon),map_pcld_Cb(klon),map_tcld_Cb(klon))
     137      allocate (map_emis_ThCi(klon),map_pcld_ThCi(klon),map_tcld_ThCi(klon))
     138      allocate (map_emis_Anv(klon),map_pcld_Anv(klon),map_tcld_Anv(klon))
     139!     endif
     140
    105141    IF (ok_hines) allocate(zustr_gwd_hines(klon), zvstr_gwd_hines(klon))
    106142    IF (.not.ok_hines.and.ok_gwd_rando) &
    … …  
    115151    IMPLICIT NONE
    116152
     153    include "clesphys.h"
     154
    117155    deallocate(snow_o,zfra_o,itau_con)
    118156    deallocate (bils_ec,bils_ech,bils_tke,bils_diss,bils_kinetic,bils_enthalp,bils_latent)
     157
     158! Marine
     159! Variables de sortie simulateur AIRS
     160
     161 !    if (ok_airs) then
     162      deallocate (map_prop_hc,map_prop_hist)
     163      deallocate (alt_tropo)
     164      deallocate (map_emis_hc,map_iwp_hc,map_deltaz_hc)
     165      deallocate (map_pcld_hc,map_tcld_hc)
     166      deallocate (map_emis_hist,map_iwp_hist,map_deltaz_hist)
     167      deallocate (map_rad_hist)
     168      deallocate (map_ntot,map_hc,map_hist)
     169      deallocate (map_Cb,map_ThCi,map_Anv)
     170      deallocate (map_emis_Cb,map_pcld_Cb,map_tcld_Cb)
     171      deallocate (map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi)
     172      deallocate (map_emis_Anv,map_pcld_Anv,map_tcld_Anv)
     173  !   endif
    119174
    120175  END SUBROUTINE phys_output_var_end

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

    r2542 r2594  
    6060         o_fsw_srf, o_wbils_srf, o_wbilo_srf, &
    6161         o_tke_srf, o_tke_max_srf,o_dltpbltke_srf, o_wstar, &
     62         o_l_mixmin,o_l_mix, &
    6263         o_cdrm, o_cdrh, o_cldl, o_cldm, o_cldh, &
    6364         o_cldt, o_JrNt, o_cldljn, o_cldmjn, &
    … …  
    116117         o_lcc3dstra, o_reffclwtop, o_ec550aer, &
    117118         o_lwcon, o_iwcon, o_temp, o_theta, &
    118          o_ovapinit, o_ovap, o_oliq, o_geop, &
     119         o_ovapinit, o_ovap, o_oliq, o_ocond, o_geop, &
    119120         o_vitu, o_vitv, o_vitw, o_pres, o_paprs, &
    120121         o_zfull, o_zhalf, o_rneb, o_rnebjn, o_rnebcon, &
    … …  
    167168         o_sens_prec_sol_oce, o_sens_prec_sol_sic, &
    168169         o_lat_prec_liq_oce, o_lat_prec_liq_sic, &
    169          o_lat_prec_sol_oce, o_lat_prec_sol_sic
     170         o_lat_prec_sol_oce, o_lat_prec_sol_sic, &
     171! Marine
     172         o_map_prop_hc, o_map_prop_hist, o_map_emis_hc, o_map_iwp_hc, &
     173         o_map_deltaz_hc, o_map_pcld_hc, o_map_tcld_hc, &
     174         o_map_emis_hist, o_map_iwp_hist, o_map_deltaz_hist, &
     175         o_map_rad_hist, &
     176         o_map_emis_Cb, o_map_pcld_Cb, o_map_tcld_Cb, &
     177         o_map_emis_ThCi, o_map_pcld_ThCi, o_map_tcld_ThCi, &
     178         o_map_emis_Anv, o_map_pcld_Anv, o_map_tcld_Anv, &
     179         o_map_ntot, o_map_hc,o_map_hist,o_map_Cb,o_map_ThCi,o_map_Anv, &
     180         o_alt_tropo
     181
    170182
    171183    USE phys_state_var_mod, only: pctsrf, paire_ter, rain_fall, snow_fall, &
    … …  
    198210    USE phys_local_var_mod, only: zxfluxlat, slp, ptstar, pt0, zxtsol, zt2m, &
    199211         t2m_min_mon, t2m_max_mon, evap, &
     212         l_mixmin,l_mix, &
    200213         zu10m, zv10m, zq2m, zustar, zxqsurf, &
    201214         rain_lsc, rain_num, snow_lsc, bils, sens, fder, &
    … …  
    237250         lcc, lcc3d, lcc3dcon, lcc3dstra, reffclwtop, &
    238251         ec550aer, flwc, fiwc, t_seri, theta, q_seri, &
    239          ql_seri, tr_seri, &
     252         ql_seri, qs_seri, tr_seri, &
    240253         zphi, u_seri, v_seri, omega, cldfra, &
    241254         rneb, rnebjn, zx_rh, d_t_dyn,  &
    … …  
    262275         zustr_gwd_hines, zvstr_gwd_hines,zustr_gwd_rando, zvstr_gwd_rando, &
    263276         zustr_gwd_front, zvstr_gwd_front,     &
    264          sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
     277         sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o, &
     278! Marine
     279         map_prop_hc, map_prop_hist, &
     280         map_emis_hc,map_iwp_hc,map_deltaz_hc,&
     281         map_pcld_hc,map_tcld_hc,&
     282         map_emis_hist,map_iwp_hist,map_deltaz_hist,&
     283         map_rad_hist,&
     284         map_ntot,map_hc,map_hist,&
     285         map_Cb,map_ThCi,map_Anv,&
     286         map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
     287         map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
     288         map_emis_Anv,map_pcld_Anv,map_tcld_Anv, &
     289         alt_tropo
     290
    265291 
    266292
    … …  
    331357
    332358    ! On calcul le nouveau tau:
    333     itau_w = itau_phy + itap + start_time * day_step_phy
     359    itau_w = itau_phy + itap
    334360    ! On le donne à iophy pour que les histwrite y aient accès:
    335361    CALL set_itau_iophy(itau_w)
    … …  
    378404       CALL histwrite_phy(o_aireTER, paire_ter)
    379405!!! Champs 2D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     406! Simulateur AIRS
     407     IF (ok_airs) then
     408      CALL histwrite_phy(o_alt_tropo,alt_tropo)
     409 
     410      CALL histwrite_phy(o_map_prop_hc,map_prop_hc)
     411      CALL histwrite_phy(o_map_prop_hist,map_prop_hist)
     412
     413      CALL histwrite_phy(o_map_emis_hc,map_emis_hc)
     414      CALL histwrite_phy(o_map_iwp_hc,map_iwp_hc)
     415      CALL histwrite_phy(o_map_deltaz_hc,map_deltaz_hc)
     416      CALL histwrite_phy(o_map_pcld_hc,map_pcld_hc)
     417      CALL histwrite_phy(o_map_tcld_hc,map_tcld_hc)
     418
     419      CALL histwrite_phy(o_map_emis_hist,map_emis_hist)
     420      CALL histwrite_phy(o_map_iwp_hist,map_iwp_hist)
     421      CALL histwrite_phy(o_map_deltaz_hist,map_deltaz_hist)
     422
     423      CALL histwrite_phy(o_map_ntot,map_ntot)
     424      CALL histwrite_phy(o_map_hc,map_hc)
     425      CALL histwrite_phy(o_map_hist,map_hist)
     426
     427      CALL histwrite_phy(o_map_Cb,map_Cb)
     428      CALL histwrite_phy(o_map_ThCi,map_ThCi)
     429      CALL histwrite_phy(o_map_Anv,map_Anv)
     430
     431      CALL histwrite_phy(o_map_emis_Cb,map_emis_Cb)
     432      CALL histwrite_phy(o_map_pcld_Cb,map_pcld_Cb)
     433      CALL histwrite_phy(o_map_tcld_Cb,map_tcld_Cb)
     434
     435      CALL histwrite_phy(o_map_emis_ThCi,map_emis_ThCi)
     436      CALL histwrite_phy(o_map_pcld_ThCi,map_pcld_ThCi)
     437      CALL histwrite_phy(o_map_tcld_ThCi,map_tcld_ThCi)
     438
     439      CALL histwrite_phy(o_map_emis_Anv,map_emis_Anv)
     440      CALL histwrite_phy(o_map_pcld_Anv,map_pcld_Anv)
     441      CALL histwrite_phy(o_map_tcld_Anv,map_tcld_Anv)
     442     ENDIF
     443
    380444       CALL histwrite_phy(o_flat, zxfluxlat)
    381445       CALL histwrite_phy(o_ptstar, ptstar)
    … …  
    619683          IF (iflag_pbl > 1) THEN
    620684             CALL histwrite_phy(o_tke_srf(nsrf),  pbl_tke(:,1:klev,nsrf))
     685             CALL histwrite_phy(o_l_mix(nsrf),  l_mix(:,1:klev,nsrf))
     686             CALL histwrite_phy(o_l_mixmin(nsrf),  l_mixmin(:,1:klev,nsrf))
    621687             CALL histwrite_phy(o_tke_max_srf(nsrf),  pbl_tke(:,1:klev,nsrf))
    622688          ENDIF
    … …  
    10551121       CALL histwrite_phy(o_ovap, q_seri)
    10561122       CALL histwrite_phy(o_oliq, ql_seri)
     1123       CALL histwrite_phy(o_ocond, ql_seri+qs_seri)
    10571124       CALL histwrite_phy(o_geop, zphi)
    10581125       CALL histwrite_phy(o_vitu, u_seri)

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

    r2570 r2594  
    10771077       !RC
    10781078       ustar(:,:)=0.
    1079        u10m(:,:)=0.
    1080        v10m(:,:)=0.
     1079!       u10m(:,:)=0.
     1080!       v10m(:,:)=0.
    10811081       rain_con(:)=0.
    10821082       snow_con(:)=0.
    … …  
    11651165
    11661166       CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
     1167!jyg<
    11671168       IF (klon_glo==1) THEN
    1168           coefh=0. ; coefm=0. ; pbl_tke=0.
    1169           coefh(:,2,:)=1.e-2 ; coefm(:,2,:)=1.e-2 ; pbl_tke(:,2,:)=1.e-2
    1170           PRINT*,'FH WARNING : lignes a supprimer'
     1169          pbl_tke(:,:,is_ave) = 0.
     1170          DO nsrf=1,nbsrf
     1171            DO k = 1,klev+1
     1172                 pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
     1173                     +pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
     1174            ENDDO
     1175          ENDDO
     1176!>jyg
    11711177       ENDIF
    11721178       !IM begin
    … …  
    14241430               klon, &
    14251431               nqtot, &
     1432               nqo, &
    14261433               pdtphys, &
    14271434               annee_ref, &
    … …  
    22942301                    wake_pe, wake_fip,  &
    22952302                    Ale_bl, Ale_bl_trig, Alp_bl, &
    2296                     Ale, Alp )
     2303                    Ale, Alp , Ale_wake, Alp_wake)
    22972304!>jyg
    22982305!
    … …  
    27762783!jyg<
    27772784!
    2778           CALL alpale_th( dtime, lmax_th, t_seri, &
     2785          CALL alpale_th( dtime, lmax_th, t_seri, cell_area,  &
    27792786                          cin, s2, n2,  &
    27802787                          ale_bl_trig, ale_bl_stat, ale_bl,  &
    2781                           alp_bl, alp_bl_stat )
     2788                          alp_bl, alp_bl_stat, &
     2789                          proba_notrig, random_notrig)
     2790
     2791          ! ------------------------------------------------------------------
     2792          ! Transport de la TKE par les panaches thermiques.
     2793          ! FH : 2010/02/01
     2794          !     if (iflag_pbl.eq.10) then
     2795          !     call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
     2796          !    s           rg,paprs,pbl_tke)
     2797          !     endif
     2798          ! -------------------------------------------------------------------
    27822799
    27832800          do i=1,klon
    … …  
    39313948#endif
    39323949    ENDIF  !ok_cosp
     3950
     3951
     3952! Marine
     3953
     3954  IF (ok_airs) then
     3955
     3956  IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/dtime)).EQ.0) THEN
     3957  write(*,*) 'je vais appeler simu_airs, ok_airs, freq_airs=', &
     3958     & ok_airs, freq_airs
     3959  call simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
     3960     & map_prop_hc,map_prop_hist,&
     3961     & map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
     3962     & map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
     3963     & map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
     3964     & map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
     3965     & map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
     3966     & map_ntot,map_hc,map_hist,&
     3967     & map_Cb,map_ThCi,map_Anv,&
     3968     & alt_tropo )
     3969  ENDIF
     3970
     3971  ENDIF  ! ok_airs
     3972
     3973
    39333974    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    39343975    !AA

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

    r2542 r2594  
    5555    data piz_strat  /0.9999998, 0.99762493/
    5656    data cg_strat   /0.73107845,0.73229635/
    57     real, dimension(nwave) :: alpha_strat_wave
     57    real, dimension(nwave_sw) :: alpha_strat_wave
    5858    data alpha_strat_wave/3.36780953,3.34667683,3.20444202,3.0293026,2.82108808/
    5959
    … …  
    153153
    154154!--total vertical aod at the 6 wavelengths
    155     DO wave=1, nwave
     155    DO wave=1, nwave_sw
    156156    DO k=1, klev
    157157    tausum_aero(:,wave,id_STRAT_phy)=tausum_aero(:,wave,id_STRAT_phy)+tau_aer_strat(:,k)*alpha_strat_wave(wave)/alpha_strat_wave(2)

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

    r2211 r2594  
    7070  ! Local
    7171  !
    72   INTEGER, PARAMETER :: las = nwave
     72  INTEGER, PARAMETER :: las = nwave_sw
    7373  LOGICAL :: soluble
    7474 

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

    r2542 r2594  
    6161!--diagnostics AOD in the SW
    6262! alpha_sw_strat_wave is *not* normalised by the 550 nm extinction coefficient
    63     real, dimension(nwave) :: alpha_sw_strat_wave
     63    real, dimension(nwave_sw) :: alpha_sw_strat_wave
    6464    data alpha_sw_strat_wave/3.708007,4.125824,4.136584,3.887478,3.507738/
    6565!
    66 !--diagnostics AOD in the LW at 10 um
    67     real :: alpha_lw_strat_wave
     66!--diagnostics AOD in the LW at 10 um (not normalised by the 550 nm ext coefficient
     67    real :: alpha_lw_strat_wave(nwave_lw)
    6868    data alpha_lw_strat_wave/0.2746812/
    6969!
    … …  
    171171
    172172!--total vertical aod at the 5 SW wavelengths
    173     DO wave=1, nwave
     173    DO wave=1, nwave_sw
    174174    DO k=1, klev
    175     tausum_aero(:,wave,id_STRAT_phy)=tausum_aero(:,wave,id_STRAT_phy)+ &
    176        tau_aer_strat(:,k)*alpha_sw_strat_wave(wave)/alpha_sw_strat_wave(2)
     175      tausum_aero(:,wave,id_STRAT_phy)=tausum_aero(:,wave,id_STRAT_phy)+ &
     176          tau_aer_strat(:,k)*alpha_sw_strat_wave(wave)/alpha_sw_strat_wave(2)
    177177    ENDDO
    178178    ENDDO
    … …  
    210210    ENDIF
    211211
     212!--total vertical aod at the 1 LW wavelength
     213    DO wave=1, nwave_lw
     214    DO k=1, klev
     215      tausum_aero(:,nwave_sw+wave,id_STRAT_phy)=tausum_aero(:,nwave_sw+wave,id_STRAT_phy)+ &
     216         tau_aer_strat(:,k)*alpha_lw_strat_wave(wave)/alpha_sw_strat_wave(2)
     217    ENDDO
     218    ENDDO
     219
    212220    DO band=1, nbands_lw_rrtm
    213221    tau_aero_lw_rrtm(:,:,2,band)  = tau_aero_lw_rrtm(:,:,2,band) + alpha_lw_abs_rrtm(band)*tau_aer_strat(:,:)

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

    r2542 r2594  
    269269!--total vertical aod at the 5 SW wavelengths
    270270!--for now use band 3 AOD into all 5 wavelengths
     271!--it is only a reasonable approximation for 550 nm (wave=2)
    271272    band=3
    272273    DO i=1, klon
    273274    DO k=1, klev
    274275      IF (stratomask(i,k).GT.0.999999) THEN
    275         DO wave=1, nwave
    276           tausum_aero(:,wave,id_STRAT_phy)=tausum_aero(:,wave,id_STRAT_phy)+tau_aer_strat(:,k,band)
     276        DO wave=1, nwave_sw
     277          tausum_aero(i,wave,id_STRAT_phy)=tausum_aero(i,wave,id_STRAT_phy)+tau_aer_strat(i,k,band)
    277278        ENDDO
    278279      ENDIF
    … …  
    308309    ENDDO
    309310
     311!--total vertical aod at 10 um
     312!--this is approximated from band 7 of RRTM
     313    band=7
     314    DO i=1, klon
     315    DO k=1, klev
     316      IF (stratomask(i,k).GT.0.999999) THEN
     317        DO wave=1, nwave_lw
     318          tausum_aero(i,nwave_sw+wave,id_STRAT_phy)=tausum_aero(i,nwave_sw+wave,id_STRAT_phy)+taulw_aer_strat(i,k,band)
     319        ENDDO
     320      ENDIF
     321    ENDDO
     322    ENDDO
     323
    310324    DO band=1, NLW
    311325    WHERE (stratomask.GT.0.999999)

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

    r2435 r2594  
    2424    USE surface_data, ONLY    : ok_veget
    2525
     26    ! See comments in each module surf_land_orchidee_xxx for compatiblity with ORCHIDEE
    2627#ifdef ORCHIDEE_NOOPENMP
     28    ! Compilation with cpp key ORCHIDEE NOOPENMP
    2729    USE surf_land_orchidee_noopenmp_mod
    2830#else
     31#if ORCHIDEE_NOZ0H
     32    ! Compilation with cpp key ORCHIDEE NOZ0H
     33    USE surf_land_orchidee_noz0h_mod
     34#else
     35    ! Compilation with default interface
    2936    USE surf_land_orchidee_mod
    3037#endif
     38#endif
     39
    3140    USE surf_land_bucket_mod
    3241    USE calcul_fluxs_mod
    … …  
    141150            evap, fluxsens, fluxlat, &             
    142151            tsol_rad, tsurf_new, alb1_new, alb2_new, &
    143             emis_new, z0m, qsurf)       
    144         z0h(1:knon)=z0m(1:knon) ! En attendant mieux
     152            emis_new, z0m, z0h, qsurf)       
    145153
    146154! 

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

    r2435 r2594  
    22MODULE surf_land_orchidee_mod
    33#ifndef ORCHIDEE_NOOPENMP
     4#ifndef ORCHIDEE_NOZ0H
    45!
    56! This module controles the interface towards the model ORCHIDEE.
    67!
    78! Compatibility with ORCHIDIEE :
    8 ! The current version can be used with ORCHIDEE/trunk from revision 2961.
    9 ! This interface can also be used with ORCHIDEE/trunk revision 1078-2960 if changing
    10 ! coszang=yrmu0 into sinang=yrmu0 at 2 places later below in this module.
     9! The current version can be used with ORCHIDEE/trunk from revision 3525.
     10! This interface is used if none of the cpp keys ORCHIDEE_NOOPENMP or ORCHIDEE_NOZ0H is set.
    1111!
    1212! Subroutines in this module : surf_land_orchidee
    … …  
    4444       evap, fluxsens, fluxlat, &             
    4545       tsol_rad, tsurf_new, alb1_new, alb2_new, &
    46        emis_new, z0_new, qsurf)
     46       emis_new, z0m_new, z0h_new, qsurf)
    4747
    4848    USE mod_surf_para
    … …  
    104104!   alb2_new     albedo in near IR interval
    105105!   emis_new     emissivite
    106 !   z0_new       surface roughness
     106!   z0m_new      surface roughness for momentum
     107!   z0h_new      surface roughness for heat
    107108!   qsurf        air moisture at surface
    108109!
    … …  
    137138    REAL, DIMENSION(klon), INTENT(OUT)        :: tsol_rad, tsurf_new
    138139    REAL, DIMENSION(klon), INTENT(OUT)        :: alb1_new, alb2_new
    139     REAL, DIMENSION(klon), INTENT(OUT)        :: emis_new, z0_new
     140    REAL, DIMENSION(klon), INTENT(OUT)        :: emis_new, z0m_new, z0h_new
    140141
    141142! Local
    … …  
    403404
    404405#ifdef CPP_VEGET
    405           CALL intersurf_main (itime+itau_phy-1, nbp_lon, nbp_lat, knon, ktindex, dtime, &
    406                lrestart_read, lrestart_write, lalo, &
    407                contfrac, neighbours, resolution, date0, &
    408                zlev,  u1_lay, v1_lay, spechum, temp_air, epot_air, ccanopy, &
     406         CALL intersurf_initialize_gathered (itime+itau_phy-1, nbp_lon, nbp_lat, knon, ktindex, dtime, &
     407               lrestart_read, lrestart_write, lalo, contfrac, neighbours, resolution, date0, &
     408               zlev,  u1_lay, v1_lay, spechum, temp_air, epot_air, &
    409409               cdrag, petA_orc, peqA_orc, petB_orc, peqB_orc, &
    410410               precip_rain, precip_snow, lwdown, swnet, swdown, ps, &
    411411               evap, fluxsens, fluxlat, coastalflow, riverflow, &
    412                tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
    413                lon_scat, lat_scat, q2m, t2m, coszang=yrmu0)
     412               tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0m_new, &
     413               lon_scat, lat_scat, q2m, t2m, z0h_new)
    414414#endif         
    415415       ENDIF
    … …  
    427427    IF (knon > 0) THEN
    428428#ifdef CPP_VEGET   
    429        CALL intersurf_main (itime+itau_phy, nbp_lon, nbp_lat, knon, ktindex, dtime,  &
     429       CALL intersurf_main_gathered (itime+itau_phy, nbp_lon, nbp_lat, knon, ktindex, dtime,  &
    430430            lrestart_read, lrestart_write, lalo, &
    431431            contfrac, neighbours, resolution, date0, &
    … …  
    434434            precip_rain(1:knon), precip_snow(1:knon), lwdown(1:knon), swnet(1:knon), swdown_vrai(1:knon), ps(1:knon), &
    435435            evap(1:knon), fluxsens(1:knon), fluxlat(1:knon), coastalflow(1:knon), riverflow(1:knon), &
    436             tsol_rad(1:knon), tsurf_new(1:knon), qsurf(1:knon), albedo_out(1:knon,:), emis_new(1:knon), z0_new(1:knon), &
    437             lon_scat, lat_scat, q2m, t2m, coszang=yrmu0(1:knon))
     436            tsol_rad(1:knon), tsurf_new(1:knon), qsurf(1:knon), albedo_out(1:knon,:), emis_new(1:knon), z0m_new(1:knon), &
     437            lon_scat, lat_scat, q2m, t2m, z0h_new(1:knon), coszang=yrmu0(1:knon))
    438438#endif       
    439439    ENDIF
    … …  
    664664!
    665665#endif
     666#endif
    666667END MODULE surf_land_orchidee_mod

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

    r2435 r2594  
    3838       knindex, rlon, rlat, yrmu0, pctsrf, &
    3939       debut, lafin, &
    40        plev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, &
     40       plev,  u1_lay, v1_lay, gustiness, temp_air, spechum, epot_air, ccanopy, &
    4141       tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
    4242       precip_rain, precip_snow, lwdown, swnet, swdown, &
    … …  
    4444       evap, fluxsens, fluxlat, &             
    4545       tsol_rad, tsurf_new, alb1_new, alb2_new, &
    46        emis_new, z0_new, qsurf)
     46       emis_new, z0_new, z0h_new, qsurf)
    4747!   
    4848! Cette routine sert d'interface entre le modele atmospherique et le
    … …  
    9595!   emis_new     emissivite
    9696!   z0_new       surface roughness
     97!   z0h_new      surface roughness, it is the same as z0_new
    9798!   qsurf        air moisture at surface
    9899!
    … …  
    121122    REAL, DIMENSION(klon), INTENT(IN)         :: yrmu0
    122123    REAL, DIMENSION(klon), INTENT(IN)         :: plev
    123     REAL, DIMENSION(klon), INTENT(IN)         :: u1_lay, v1_lay
     124    REAL, DIMENSION(klon), INTENT(IN)         :: u1_lay, v1_lay, gustiness
    124125    REAL, DIMENSION(klon), INTENT(IN)         :: temp_air, spechum
    125126    REAL, DIMENSION(klon), INTENT(IN)         :: epot_air, ccanopy
    … …  
    136137    REAL, DIMENSION(klon), INTENT(OUT)        :: tsol_rad, tsurf_new
    137138    REAL, DIMENSION(klon), INTENT(OUT)        :: alb1_new, alb2_new
    138     REAL, DIMENSION(klon), INTENT(OUT)        :: emis_new, z0_new
     139    REAL, DIMENSION(klon), INTENT(OUT)        :: emis_new, z0_new, z0h_new
    139140
    140141! Local
    … …  
    496497
    497498    albedo_keep(1:knon) = (albedo_out(1:knon,1)+albedo_out(1:knon,2))/2.
     499   
     500    ! ORCHIDEE only gives one value for z0_new. Copy it into z0h_new.
     501    z0h_new(:)=z0_new(:)
    498502
    499503!* Send to coupler

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

    r2488 r2594  
    8383  USE print_control_mod, ONLY: lunout
    8484  IMPLICIT NONE
     85  INCLUDE 'YOMCST.h'
    8586    REAL,INTENT(IN) :: pdtphys_
    8687    REAL            :: julian_date
    … …  
    9495   
    9596    ! Update elapsed time since begining of run:
    96     current_time=current_time+pdtphys
     97    current_time=current_time+pdtphys/rday
    9798
    9899    ! Compute corresponding Julian date and update calendar
    99     CALL ymds2ju(annee_ref,1,day_ini,start_time+current_time,julian_date)
     100    CALL ymds2ju(annee_ref,1,day_ini,(start_time+current_time)*rday,julian_date)
    100101    CALL phys_cal_update(julian_date)
    101102   

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

    r2471 r2594  
    1 
    2 ! $Header$
    3 
    4 SUBROUTINE yamada4(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
     1!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2
     3SUBROUTINE yamada4(ni, nsrf, ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
    54    cd, q2, km, kn, kq, ustar, iflag_pbl)
     5
    66  USE dimphy
    7   USE print_control_mod, ONLY: prt_level
    87  USE ioipsl_getin_p_mod, ONLY : getin_p
    98
    109  IMPLICIT NONE
    11 
     10  include "iniprint.h"
     11  ! .......................................................................
     12  ! ym#include "dimensions.h"
     13  ! ym#include "dimphy.h"
     14  ! ************************************************************************************************
     15  !
     16  ! yamada4: subroutine qui calcule le transfert turbulent avec une fermeture d'ordre 2 ou 2.5
     17  !
     18  ! Reference: Simulation of nocturnal drainage flows by a q2l Turbulence Closure Model
     19  !            Yamada T.
     20  !            J. Atmos. Sci, 40, 91-106, 1983
     21  !
     22  !************************************************************************************************
     23  ! Input :
     24  !'======
     25  ! ni: indice horizontal sur la grille de base, non restreinte
     26  ! nsrf: type de surface
     27  ! ngrid: nombre de mailles concern??es sur l'horizontal
    1228  ! dt : pas de temps
    1329  ! g  : g
     30  ! rconst: constante de l'air sec
    1431  ! zlev : altitude a chaque niveau (interface inferieure de la couche
    1532  ! de meme indice)
    … …  
    1734  ! u,v : vitesse au centre de chaque couche
    1835  ! (en entree : la valeur au debut du pas de temps)
    19   ! teta : temperature potentielle au centre de chaque couche
     36  ! teta : temperature potentielle virtuelle au centre de chaque couche
    2037  ! (en entree : la valeur au debut du pas de temps)
    21   ! cd : cdrag
     38  ! cd : cdrag pour la quantit?? de mouvement
    2239  ! (en entree : la valeur au debut du pas de temps)
     40  ! ustar: vitesse de friction calcul??e par une formule diagnostique
     41  ! iflag_pbl: flag pour choisir des options du sch??ma de turbulence
     42  !
     43  !             iflag_pbl doit valoir entre 6 et 9
     44  !             l=6, on prend  systematiquement une longueur d'equilibre
     45  !             iflag_pbl=6 : MY 2.0
     46  !             iflag_pbl=7 : MY 2.0.Fournier
     47  !             iflag_pbl=8/9 : MY 2.5
     48  !             iflag_pbl=8 with special obsolete treatments for convergence
     49  !             with Cmpi5 NPv3.1 simulations
     50  !             iflag_pbl=10/11 :  New scheme M2 and N2 explicit and dissiptation exact
     51  !             iflag_pbl=12 = 11 with vertical diffusion off q2
     52  !
     53  !             2013/04/01 (FH hourdin@lmd.jussieu.fr)
     54  !             Correction for very stable PBLs (iflag_pbl=10 and 11)
     55  !             iflag_pbl=8 converges numerically with NPv3.1
     56  !             iflag_pbl=11 -> the model starts with NP from start files created by ce0l
     57  !                          -> the model can run with longer time-steps.
     58  !
     59  ! Inpout/Output :
     60  !==============
    2361  ! q2 : $q^2$ au bas de chaque couche
    2462  ! (en entree : la valeur au debut du pas de temps)
    2563  ! (en sortie : la valeur a la fin du pas de temps)
     64 
     65  ! Outputs:
     66  !==========
    2667  ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
    2768  ! couche)
    … …  
    2970  ! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
    3071  ! (en sortie : la valeur a la fin du pas de temps)
    31 
    32   ! iflag_pbl doit valoir entre 6 et 9
    33   ! l=6, on prend  systematiquement une longueur d'equilibre
    34   ! iflag_pbl=6 : MY 2.0
    35   ! iflag_pbl=7 : MY 2.0.Fournier
    36   ! iflag_pbl=8/9 : MY 2.5
    37   ! iflag_pbl=8 with special obsolete treatments for convergence
    38   ! with Cmpi5 NPv3.1 simulations
    39   ! iflag_pbl=10/11 :  New scheme M2 and N2 explicit and dissiptation exact
    40   ! iflag_pbl=12 = 11 with vertical diffusion off q2
    41 
    42   ! 2013/04/01 (FH hourdin@lmd.jussieu.fr)
    43   ! Correction for very stable PBLs (iflag_pbl=10 and 11)
    44   ! iflag_pbl=8 converges numerically with NPv3.1
    45   ! iflag_pbl=11 -> the model starts with NP from start files created by ce0l
    46   ! -> the model can run with longer time-steps.
    47   ! .......................................................................
    48 
     72  !
     73  !.......................................................................
     74
     75  !=======================================================================
     76  ! Declarations:
     77  !=======================================================================
     78
     79
     80  ! Inputs/Outputs
     81  !----------------
    4982  REAL dt, g, rconst
    5083  REAL plev(klon, klev+1), temp(klon, klev)
    … …  
    5790  REAL teta(klon, klev)
    5891  REAL cd(klon)
    59   REAL q2(klon, klev+1), qpre
     92  REAL q2(klon, klev+1)
    6093  REAL unsdz(klon, klev)
    6194  REAL unsdzdec(klon, klev+1)
    62 
     95  REAL kn(klon, klev+1)
    6396  REAL km(klon, klev+1)
    64   REAL kmpre(klon, klev+1), tmp2
     97  INTEGER iflag_pbl, ngrid, nsrf
     98  INTEGER ni(klon)
     99
     100
     101  ! Local
     102  !-------
     103
     104  INCLUDE "clesphys.h"
     105
     106
     107  REAL kmpre(klon, klev+1), tmp2, qpre
    65108  REAL mpre(klon, klev+1)
    66   REAL kn(klon, klev+1)
    67109  REAL kq(klon, klev+1)
    68110  REAL ff(klon, klev+1), delta(klon, klev+1)
    69111  REAL aa(klon, klev+1), aa0, aa1
    70   INTEGER iflag_pbl, ngrid
    71112  INTEGER nlay, nlev
    72 
    73113  LOGICAL first
    74114  INTEGER ipas
    … …  
    77117  DATA first, ipas/.FALSE., 0/
    78118  !$OMP THREADPRIVATE( first,ipas)
    79   REAL,SAVE :: lmixmin=1.
    80   !$OMP THREADPRIVATE(lmixmin)
    81 
    82 
     119  LOGICAL hboville
    83120  INTEGER ig, k
    84 
    85 
    86121  REAL ri, zrif, zalpha, zsm, zsn
    87122  REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
    88 
    89123  REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
    90124  REAL dtetadz(klon, klev+1)
    91125  REAL m2cstat, mcstat, kmcstat
    92126  REAL l(klon, klev+1)
    93   REAL, ALLOCATABLE, SAVE :: l0(:)
    94   !$OMP THREADPRIVATE(l0)
    95   REAL sq(klon), sqz(klon), zz(klon, klev+1)
     127  REAL zz(klon, klev+1)
    96128  INTEGER iter
    97 
    98129  REAL ric, rifc, b1, kap
    99130  SAVE ric, rifc, b1, kap
    100131  DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/
    101132  !$OMP THREADPRIVATE(ric,rifc,b1,kap)
     133  REAL seuilsm, seuilalpha
     134  REAL,SAVE :: lmixmin
     135  !$OMP THREADPRIVATE(lmixmin)
    102136  REAL frif, falpha, fsm
    103   REAL fl, zzz, zl0, zq2, zn2
    104 
    105137  REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), &
    106138    lyam(klon, klev), knyam(klon, klev), w2yam(klon, klev), t2yam(klon, klev)
    107139  LOGICAL, SAVE :: firstcall = .TRUE.
    108140  !$OMP THREADPRIVATE(firstcall)
     141
     142
     143  ! Fonctions utilis??es
     144  !--------------------
     145
    109146  frif(ri) = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
    110147  falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri)
    111148  fsm(ri) = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
    112   fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)*kap*zlev(ig, &
    113     k)/(kap*zlev(ig,k)+l0(ig)),0.5*sqrt(q2(ig,k))/sqrt( &
    114     max(n2(ig,k),1.E-10))), lmixmin)
    115 
    116 
    117   nlay = klev
    118   nlev = klev + 1
     149 
    119150
    120151  IF (firstcall) THEN
    121     ALLOCATE (l0(klon))
    122152    firstcall = .FALSE.
     153    lmixmin=1.
    123154    CALL getin_p('lmixmin',lmixmin)
    124155  END IF
     156
     157
     158
     159!===============================================================================
     160! Flags, tests et ??valuations de constantes
     161!===============================================================================
     162
     163! On utilise ou non la routine de Holstalg Boville pour les cas tres stables
     164   hboville=.TRUE.
    125165
    126166
    … …  
    129169  END IF
    130170
     171! Seuil dans le code de turbulence
     172 seuilalpha=1.12
     173 seuilsm=0.085
     174
     175  nlay = klev
     176  nlev = klev + 1
    131177  ipas = ipas + 1
    132178
    133179
    134   ! .......................................................................
    135   ! les increments verticaux
    136   ! .......................................................................
    137 
    138   ! !!!!! allerte !!!!!c
    139   ! !!!!! zlev n'est pas declare a nlev !!!!!c
    140   ! !!!!! ---->
     180!========================================================================
     181! Calcul des increments verticaux
     182!=========================================================================
     183
     184 
     185! Attention: zlev n'est pas declare a nlev
    141186  DO ig = 1, ngrid
    142187    zlev(ig, nlev) = zlay(ig, nlay) + (zlay(ig,nlay)-zlev(ig,nlev-1))
    143188  END DO
    144   ! !!!!! <----
    145   ! !!!!! allerte !!!!!c
     189
    146190
    147191  DO k = 1, nlay
    … …  
    162206  END DO
    163207
    164   ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    165   ! Computing M^2, N^2, Richardson numbers, stability functions
    166   ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    167     ! initialize arrays:
    168   m2(:, :) = 0.0
    169   sm(:, :) = 0.0
    170   rif(:, :) = 0.0
    171 
     208!=========================================================================
     209! Richardson number and stability functions
     210!=========================================================================
     211 
     212! initialize arrays:
     213
     214  m2(1:ngrid, :) = 0.0
     215  sm(1:ngrid, :) = 0.0
     216  rif(1:ngrid, :) = 0.0
     217
     218!------------------------------------------------------------
    172219  DO k = 2, klev
     220
    173221    DO ig = 1, ngrid
    174222      dz(ig, k) = zlay(ig, k) - zlay(ig, k-1)
    … …  
    177225      dtetadz(ig, k) = (teta(ig,k)-teta(ig,k-1))/dz(ig, k)
    178226      n2(ig, k) = g*2.*dtetadz(ig, k)/(teta(ig,k-1)+teta(ig,k))
    179       ! n2(ig,k)=0.
    180227      ri = n2(ig, k)/max(m2(ig,k), 1.E-10)
    181228      IF (ri<ric) THEN
    … …  
    184231        rif(ig, k) = rifc
    185232      END IF
     233
    186234      IF (rif(ig,k)<0.16) THEN
    187235        alpha(ig, k) = falpha(rif(ig,k))
    188236        sm(ig, k) = fsm(rif(ig,k))
    189237      ELSE
    190         alpha(ig, k) = 1.12
    191         sm(ig, k) = 0.085
     238        alpha(ig, k) = seuilalpha
     239        sm(ig, k) = seuilsm
    192240      END IF
    193241      zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k)
    … …  
    196244
    197245
    198   ! ====================================================================
    199   ! Computing the mixing length
    200   ! ====================================================================
    201 
    202   ! Mise a jour de l0
    203   IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
    204 
    205     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    206     ! Iterative computation of l0
    207     ! This version is kept for iflag_pbl only for convergence
    208     ! with NPv3.1 Cmip5 simulations
    209     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    210 
    211     DO ig = 1, ngrid
    212       sq(ig) = 1.E-10
    213       sqz(ig) = 1.E-10
    214     END DO
    215     DO k = 2, klev - 1
    216       DO ig = 1, ngrid
    217         zq = sqrt(q2(ig,k))
    218         sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1))
    219         sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
    220       END DO
    221     END DO
    222     DO ig = 1, ngrid
    223       l0(ig) = 0.2*sqz(ig)/sq(ig)
    224     END DO
     246
     247! ====================================================================
     248! Computing the mixing length
     249! ====================================================================
     250
     251 
     252  CALL mixinglength(ni,nsrf,ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, l)
     253
     254
     255
     256  !=======================================================================
     257  !     DIFFERENT TYPES  DE SCHEMA de  YAMADA
     258  !=======================================================================
     259
     260  !--------------
     261  ! Yamada 2.0
     262  !--------------
     263  IF (iflag_pbl==6) THEN
     264 
     265    DO k = 2, klev
     266      q2(1:ngrid, k) = l(1:ngrid, k)**2*zz(1:ngrid, k)
     267    END DO
     268
     269
     270  !------------------
     271  ! Yamada 2.Fournier
     272  !------------------
     273
     274  ELSE IF (iflag_pbl==7) THEN
     275
     276
     277    ! Calcul de l,  km, au pas precedent
     278    !....................................
    225279    DO k = 2, klev
    226280      DO ig = 1, ngrid
    227         l(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
    228       END DO
    229     END DO
    230     ! print*,'L0 cas 8 ou 10 ',l0
    231 
    232   ELSE
    233 
    234     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    235     ! In all other case, the assymptotic mixing length l0 is imposed (100m)
    236     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    237 
    238     l0(:) = 150.
    239     DO k = 2, klev
    240       DO ig = 1, ngrid
    241         l(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
    242       END DO
    243     END DO
    244     ! print*,'L0 cas autres ',l0
    245 
    246   END IF
    247 
    248 
    249   ! ====================================================================
    250   ! Yamada 2.0
    251   ! ====================================================================
    252   IF (iflag_pbl==6) THEN
    253 
    254     DO k = 2, klev
    255       q2(:, k) = l(:, k)**2*zz(:, k)
    256     END DO
    257 
    258 
    259   ELSE IF (iflag_pbl==7) THEN
    260     ! ====================================================================
    261     ! Yamada 2.Fournier
    262     ! ====================================================================
    263 
    264     ! Calcul de l,  km, au pas precedent
    265     DO k = 2, klev
    266       DO ig = 1, ngrid
    267         ! print*,'SMML=',sm(ig,k),l(ig,k)
    268281        delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k))
    269282        kmpre(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
    270283        mpre(ig, k) = sqrt(m2(ig,k))
    271         ! print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
    272284      END DO
    273285    END DO
    … …  
    278290        mcstat = sqrt(m2cstat)
    279291
    280         ! print*,'M2 L=',k,mpre(ig,k),mcstat
    281 
    282         ! -----{puis on ecrit la valeur de q qui annule l'equation de m
    283         ! supposee en q3}
     292     ! Puis on ecrit la valeur de q qui annule l'equation de m supposee en q3
     293     !.........................................................................
    284294
    285295        IF (k==2) THEN
    … …  
    293303            (unsdz(ig,k)+unsdz(ig,k-1))
    294304        END IF
    295         ! print*,'T2 L=',k,tmp2
     305
    296306        tmp2 = kmcstat/(sm(ig,k)/q2(ig,k))/l(ig, k)
    297307        q2(ig, k) = max(tmp2, 1.E-12)**(2./3.)
    298         ! print*,'Q2 L=',k,q2(ig,k)
    299308
    300309      END DO
    301310    END DO
    302311
     312
     313    ! ------------------------
     314    ! Yamada 2.5 a la Didi
     315    !-------------------------
     316
    303317  ELSE IF (iflag_pbl==8 .OR. iflag_pbl==9) THEN
    304     ! ====================================================================
    305     ! Yamada 2.5 a la Didi
    306     ! ====================================================================
    307 
    308 
    309     ! Calcul de l,  km, au pas precedent
     318
     319    ! Calcul de l, km, au pas precedent
     320    !....................................
    310321    DO k = 2, klev
    311322      DO ig = 1, ngrid
    312         ! print*,'SMML=',sm(ig,k),l(ig,k)
    313323        delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k))
    314324        IF (delta(ig,k)<1.E-20) THEN
    315           ! print*,'ATTENTION   L=',k,'   Delta=',delta(ig,k)
    316325          delta(ig, k) = 1.E-20
    317326        END IF
    … …  
    319328        aa0 = (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
    320329        aa1 = (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
    321         ! abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
    322330        aa(ig, k) = aa1*dt/(delta(ig,k)*l(ig,k))
    323         ! print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
    324331        qpre = sqrt(q2(ig,k))
    325         ! if (iflag_pbl.eq.8 ) then
    326332        IF (aa(ig,k)>0.) THEN
    327333          q2(ig, k) = (qpre+aa(ig,k)*qpre*qpre)**2
    … …  
    337343        ! endif
    338344        q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
    339         ! print*,'Q2 L=',k,q2(ig,k),qpre*qpre
    340345      END DO
    341346    END DO
    … …  
    343348  ELSE IF (iflag_pbl>=10) THEN
    344349
    345     ! print*,'Schema mixte D'
    346     ! print*,'Longueur ',l(:,:)
    347350    DO k = 2, klev - 1
    348       DO ig = 1, ngrid
    349         l(ig, k) = max(l(ig,k), lmixmin)
    350         km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
    351         q2(ig, k) = q2(ig, k) + dt*km(ig, k)*m2(ig, k)*(1.-rif(ig,k))
    352         q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
    353         q2(ig, k) = 1./(1./sqrt(q2(ig,k))+dt/(2*l(ig,k)*b1))
    354         q2(ig, k) = q2(ig, k)*q2(ig, k)
    355       END DO
     351      km(1:ngrid, k) = l(1:ngrid, k)*sqrt(q2(1:ngrid,k))*sm(1:ngrid, k)
     352      q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
     353      q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4)
     354      q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1))
     355      q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k)
    356356    END DO
    357357
    … …  
    362362  END IF ! Fin du cas 8
    363363
    364   ! print*,'OK8'
    365364
    366365  ! ====================================================================
    367   ! Calcul des coefficients de m�ange
     366  ! Calcul des coefficients de melange
    368367  ! ====================================================================
     368
    369369  DO k = 2, klev
    370     ! print*,'k=',k
    371370    DO ig = 1, ngrid
    372       ! abde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
    373371      zq = sqrt(q2(ig,k))
    374       km(ig, k) = l(ig, k)*zq*sm(ig, k)
    375       kn(ig, k) = km(ig, k)*alpha(ig, k)
    376       kq(ig, k) = l(ig, k)*zq*0.2
    377       ! print*,'KML=',km(ig,k),kn(ig,k)
    378     END DO
    379   END DO
     372      km(ig, k) = l(ig, k)*zq*sm(ig, k)     ! For momentum
     373      kn(ig, k) = km(ig, k)*alpha(ig, k)    ! For scalars
     374      kq(ig, k) = l(ig, k)*zq*0.2           ! For TKE
     375    END DO
     376  END DO
     377
     378
     379  !====================================================================
     380  ! Transport diffusif vertical de la TKE par la TKE
     381  !====================================================================
     382
     383
    380384    ! initialize near-surface and top-layer mixing coefficients
    381   kq(1:ngrid, 1) = kq(1:ngrid, 2) ! constant (ie no gradient) near the surface
    382   kq(1:ngrid, klev+1) = 0 ! zero at the top
    383 
    384   ! Transport diffusif vertical de la TKE.
     385    !...........................................................
     386
     387  kq(1:ngrid, 1) = kq(1:ngrid, 2)    ! constant (ie no gradient) near the surface
     388  kq(1:ngrid, klev+1) = 0            ! zero at the top
     389
     390    ! Transport diffusif vertical de la TKE.
     391    !.......................................
     392
    385393  IF (iflag_pbl>=12) THEN
    386     ! print*,'YAMADA VDIF'
    387     q2(:, 1) = q2(:, 2)
     394    q2(1:ngrid, 1) = q2(1:ngrid, 2)
    388395    CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2)
    389396  END IF
    390397
    391   ! Traitement des cas noctrunes avec l'introduction d'une longueur
    392   ! minilale.
    393 
    394   ! ====================================================================
    395   ! Traitement particulier pour les cas tres stables.
    396   ! D'apres Holtslag Boville.
     398
     399  !====================================================================
     400  ! Traitement particulier pour les cas tres stables, introduction d'une
     401  ! longueur de m??lange minimale
     402  !====================================================================
     403  !
     404  ! Reference: Local versus Nonlocal boundary-layer diffusion in a global climate model
     405  !            Holtslag A.A.M. and Boville B.A.
     406  !            J. Clim., 6, 1825-1842, 1993
     407
     408
     409 IF (hboville) THEN
     410
    397411
    398412  IF (prt_level>1) THEN
    399413    PRINT *, 'YAMADA4 0'
    400   END IF !(prt_level>1) THEN
     414  END IF
     415
    401416  DO ig = 1, ngrid
    402417    coriol(ig) = 1.E-4
    … …  
    404419  END DO
    405420
    406   ! print*,'pblhmin ',pblhmin
    407   ! Test a remettre 21 11 02
    408   ! test abd 13 05 02      if(0.eq.1) then
    409421  IF (1==1) THEN
    410422    IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
    … …  
    419431          END IF
    420432          IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
    421             ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
    422             ! s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
     433
    423434            kn(ig, k) = kmin
    424435            km(ig, k) = kmin
    425436            kq(ig, k) = kmin
    426             ! la longueur de melange est suposee etre l= kap z
    427             ! K=l q Sm d'ou q2=(K/l Sm)**2
     437
     438 ! la longueur de melange est suposee etre l= kap z
     439 ! K=l q Sm d'ou q2=(K/l Sm)**2
     440
    428441            q2(ig, k) = (qmin/sm(ig,k))**2
    429442          END IF
    … …  
    432445
    433446    ELSE
    434 
    435447      DO k = 2, klev
    436448        DO ig = 1, ngrid
    … …  
    442454          END IF
    443455          IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
    444             ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
    445             ! s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
    446456            kn(ig, k) = kmin
    447457            km(ig, k) = kmin
    448458            kq(ig, k) = kmin
    449             ! la longueur de melange est suposee etre l= kap z
    450             ! K=l q Sm d'ou q2=(K/l Sm)**2
     459 ! la longueur de melange est suposee etre l= kap z
     460 ! K=l q Sm d'ou q2=(K/l Sm)**2
    451461            sm(ig, k) = 1.
    452462            alpha(ig, k) = 1.
    … …  
    463473  END IF
    464474
     475 END IF ! hboville
     476
    465477  IF (prt_level>1) THEN
    466478    PRINT *, 'YAMADA4 1'
    467479  END IF !(prt_level>1) THEN
     480
     481
     482 !======================================================
     483 ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane
     484 !======================================================
     485 !
     486 ! Reference: A New Second-Order Turbulence Closure Scheme for the Planetary Boundary Layer
     487 !            Abdella K and McFarlane N
     488 !            J. Atmos. Sci., 54, 1850-1867, 1997
     489
    468490  ! Diagnostique pour stokage
     491  !..........................
    469492
    470493  IF (1==0) THEN
    … …  
    482505    knyam(1:ngrid, 2:klev) = kn(1:ngrid, 2:klev)
    483506
    484     ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane
     507
     508  ! Calcul de w'2 et T'2
     509  !.......................
    485510
    486511    w2yam(1:ngrid, 2:klev) = q2(1:ngrid, 2:klev)*0.24 + &
    … …  
    493518  END IF
    494519
    495   ! print*,'OKFIN'
     520!============================================================================
     521
    496522  first = .FALSE.
    497523  RETURN
     524
     525
    498526END SUBROUTINE yamada4
     527
     528!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     529
     530!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    499531SUBROUTINE vdif_q2(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
     532
    500533  USE dimphy
    501534  IMPLICIT NONE
    502 
    503   ! dt : pas de temps
     535 
     536  include "dimensions.h"
     537
     538!    vdif_q2: subroutine qui calcule la diffusion de la TKE par la TKE
     539!             avec un schema implicite en temps avec
     540!             inversion d'un syst??me tridiagonal
     541!
     542!     Reference: Description of the interface with the surface and
     543!                the computation of the turbulet diffusion in LMDZ
     544!                Technical note on LMDZ
     545!                Dufresne, J-L, Ghattas, J. and Grandpeix, J-Y
     546!
     547!============================================================================
     548! Declarations
     549!============================================================================
    504550
    505551  REAL plev(klon, klev+1)
    … …  
    516562  INTEGER i, k
    517563
    518   ! print*,'RD=',rconst
     564
     565!=========================================================================
     566! Calcul
     567!=========================================================================
     568
    519569  DO k = 1, klev
    520570    DO i = 1, ngrid
    521       ! test
    522       ! print*,'i,k',i,k
    523       ! print*,'temp(i,k)=',temp(i,k)
    524       ! print*,'(plev(i,k)-plev(i,k+1))=',plev(i,k),plev(i,k+1)
    525571      zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
    526572      kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ &
    … …  
    546592      denom(i, k) = deltap(i, k) + (1.-beta(i,k+1))*kstar(i, k) + &
    547593        kstar(i, k-1)
    548       ! correction d'un bug 10 01 2001
    549594      alpha(i, k) = (q2(i,k)*deltap(i,k)+kstar(i,k)*alpha(i,k+1))/denom(i, k)
    550595      beta(i, k) = kstar(i, k-1)/denom(i, k)
    … …  
    553598
    554599  ! Si on recalcule q2(1)
     600  !.......................
    555601  IF (1==0) THEN
    556602    DO i = 1, ngrid
    … …  
    559605    END DO
    560606  END IF
    561   ! sinon, on peut sauter cette boucle...
     607
    562608
    563609  DO k = 2, klev + 1
    … …  
    569615  RETURN
    570616END SUBROUTINE vdif_q2
    571 SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
    572   USE dimphy
     617!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     618
     619
     620
     621!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     622 SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
     623 
     624   USE dimphy
    573625  IMPLICIT NONE
    574626
    575   ! dt : pas de temps
     627  include "dimensions.h"
     628!
     629! vdif_q2e: subroutine qui calcule la diffusion de TKE par la TKE
     630!           avec un schema explicite en temps
     631
     632
     633!====================================================
     634! Declarations
     635!====================================================
    576636
    577637  REAL plev(klon, klev+1)
    … …  
    585645  REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
    586646  INTEGER ngrid
    587 
    588647  INTEGER i, k
     648
     649
     650!==================================================
     651! Calcul
     652!==================================================
    589653
    590654  DO k = 1, klev
    … …  
    621685  RETURN
    622686END SUBROUTINE vdif_q2e
     687
     688!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     689
     690
     691!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     692
     693SUBROUTINE mixinglength(ni, nsrf, ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, lmix)
     694
     695
     696
     697  USE dimphy
     698  USE phys_state_var_mod, only: zstd, zsig, zmea
     699  USE phys_local_var_mod, only: l_mixmin, l_mix
     700
     701 ! zstd: ecart type de la'altitud e sous-maille
     702 ! zmea: altitude moyenne sous maille
     703 ! zsig: pente moyenne de le maille
     704
     705  USE geometry_mod, only: cell_area
     706  ! aire_cell: aire de la maille
     707
     708  IMPLICIT NONE
     709!*************************************************************************
     710! Subrourine qui calcule la longueur de m??lange dans le sch??ma de turbulence
     711! avec la formule de Blackadar
     712! Calcul d'un  minimum en fonction de l'orographie sous-maille:
     713! L'id??e est la suivante: plus il y a de relief, plus il y a du m??lange
     714! induit par les circulations meso et submeso ??chelles.
     715!
     716! References: * The vertical distribution of wind and turbulent exchange in a neutral atmosphere
     717!               Blackadar A.K.
     718!               J. Geophys. Res., 64, No 8, 1962
     719!
     720!             * An evaluation of neutral and convective planetary boundary-layer parametrisations relative
     721!               to large eddy simulations
     722!               Ayotte K et al
     723!               Boundary Layer Meteorology, 79, 131-175, 1996
     724!
     725!
     726!             * Local Similarity in the Stable Boundary Layer and Mixing length Approaches: consistency of concepts
     727!               Van de Wiel B.J.H et al
     728!               Boundary-Lay Meteorol, 128, 103-166, 2008
     729!
     730!
     731! Histoire:
     732!----------
     733! * premi??re r??daction, Etienne et Frederic, 09/06/2016
     734!
     735! ***********************************************************************
     736
     737!==================================================================
     738! Declarations
     739!==================================================================
     740
     741! Inputs
     742!-------
     743 INTEGER            ni(klon)           ! indice sur la grille original (non restreinte)
     744 INTEGER            nsrf               ! Type de surface
     745 INTEGER            ngrid              ! Nombre de points concern??s sur l'horizontal
     746 INTEGER            iflag_pbl          ! Choix du sch??ma de turbulence
     747 REAL            pbl_lmixmin_alpha  ! on active ou non le calcul de la longueur de melange minimum
     748 REAL               lmixmin            ! Minimum absolu de la longueur de m??lange
     749 REAL               zlay(klon, klev)   ! altitude du centre de la couche
     750 REAL               zlev(klon, klev+1) ! atitude de l'interface inf??rieure de la couche
     751 REAL               u(klon, klev)      ! vitesse du vent zonal
     752 REAL               v(klon, klev)      ! vitesse du vent meridional
     753 REAL               q2(klon, klev+1)   ! energie cin??tique turbulente
     754 REAL               n2(klon, klev+1)   ! frequence de Brunt-Vaisala
     755
     756!In/out
     757!-------
     758
     759  LOGICAL, SAVE :: firstcall = .TRUE.
     760  !$OMP THREADPRIVATE(firstcall)
     761
     762! Outputs
     763!---------
     764
     765 REAL               lmix(klon, klev+1)    ! Longueur de melange 
     766
     767
     768! Local
     769!-------
     770 
     771 INTEGER  ig,jg, k
     772 REAL     h_oro(klon)
     773 REAL     hlim(klon)
     774 REAL, SAVE :: kap=0.4,kapb=0.4
     775 REAL zq
     776 REAL sq(klon), sqz(klon)
     777 REAL, ALLOCATABLE, SAVE :: l0(:)
     778  !$OMP THREADPRIVATE(l0)
     779 REAL fl, zzz, zl0, zq2, zn2
     780 REAL famorti, zzzz, zh_oro, zhlim
     781 REAL l1(klon, klev+1), l2(klon,klev+1)
     782 REAL winds(klon, klev)
     783 REAL xcell
     784 REAL zstdslope(klon) 
     785 REAL lmax
     786 REAL l2strat, l2neutre, extent 
     787 REAL l2limit(klon)
     788!===============================================================
     789! Fonctions utiles
     790!===============================================================
     791
     792! Calcul de l suivant la formule de Blackadar 1962 adapt??e par Ayotte 1996
     793!..........................................................................
     794
     795 fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)*kap*zlev(ig, &
     796    k)/(kap*zlev(ig,k)+l0(ig)),0.5*sqrt(q2(ig,k))/sqrt( &
     797    max(n2(ig,k),1.E-10))), 1.E-5)
     798 
     799! Fonction d'amortissement de la turbulence au dessus de la montagne
     800! On augmente l'amortissement en diminuant la valeur de hlim (extent) dans le code
     801!.....................................................................
     802
     803 famorti(zzzz, zh_oro, zhlim)=(-1.)*ATAN((zzzz-zh_oro)/(zhlim-zh_oro))*2./3.1416+1.   
     804
     805  IF (ngrid==0) RETURN
     806
     807  IF (firstcall) THEN
     808    ALLOCATE (l0(klon))
     809    firstcall = .FALSE.
     810  END IF
     811
     812
     813!=====================================================================
     814!         CALCUL de la LONGUEUR de m??lange suivant BLACKADAR: l1
     815!=====================================================================
     816
     817
     818  IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
     819
     820   
     821    ! Iterative computation of l0
     822    ! This version is kept for iflag_pbl only for convergence
     823    ! with NPv3.1 Cmip5 simulations
     824    !...................................................................
     825
     826    DO ig = 1, ngrid
     827      sq(ig) = 1.E-10
     828      sqz(ig) = 1.E-10
     829    END DO
     830    DO k = 2, klev - 1
     831      DO ig = 1, ngrid
     832        zq = sqrt(q2(ig,k))
     833        sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1))
     834        sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
     835      END DO
     836    END DO
     837    DO ig = 1, ngrid
     838      l0(ig) = 0.2*sqz(ig)/sq(ig)
     839    END DO
     840    DO k = 2, klev
     841      DO ig = 1, ngrid
     842        l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
     843      END DO
     844    END DO
     845
     846  ELSE
     847
     848   
     849    ! In all other case, the assymptotic mixing length l0 is imposed (150m)
     850    !......................................................................
     851
     852    l0(1:ngrid) = 150.
     853    DO k = 2, klev
     854      DO ig = 1, ngrid
     855        l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
     856      END DO
     857    END DO
     858
     859  END IF
     860
     861!=================================================================================
     862!  CALCUL d'une longueur de melange en fonctions de la topographie sous maille: l2
     863! si plb_lmixmin_alpha=TRUE et si on se trouve sur de la terre ( pas actif sur les
     864! glacier, la glace de mer et les oc??ans)
     865!=================================================================================
     866
     867   l2(1:ngrid,:)=0.0
     868   l_mixmin(1:ngrid,:,nsrf)=0.
     869   l_mix(1:ngrid,:,nsrf)=0.
     870
     871   IF (nsrf .EQ. 1) THEN
     872
     873! coefficients
     874!--------------
     875
     876     extent=2.                                                         ! On ??tend l'impact du relief jusqu'?? extent*h, extent >1. 
     877     lmax=150.                                                         ! Longueur de m??lange max dans l'absolu
     878
     879! calculs
     880!---------
     881
     882     DO ig=1,ngrid
     883
     884      ! On calcule la hauteur du relief
     885      !.................................
     886      ! On ne peut pas prendre zstd seulement pour caracteriser le relief sous maille
     887      ! car sur un terrain pentu mais sans relief, zstd est non nul (comme en Antarctique, C. Genthon)
     888      ! On corrige donc zstd avec l'ecart type de l'altitude dans la maille sans relief
     889      ! (en gros, une maille de taille xcell avec une pente constante zstdslope)
     890      jg=ni(ig)
     891!     IF (zsig(jg) .EQ. 0.) THEN
     892!          zstdslope(ig)=0.         
     893!     ELSE
     894!     xcell=sqrt(cell_area(jg))
     895!     zstdslope(ig)=max((xcell*zsig(jg)-zmea(jg))**3 /(3.*zsig(jg)),0.)
     896!     zstdslope(ig)=sqrt(zstdslope(ig))
     897!     END IF
     898     
     899!     h_oro(ig)=max(zstd(jg)-zstdslope(ig),0.)   ! Hauteur du relief
     900      h_oro(ig)=zstd(jg)
     901      hlim(ig)=extent*h_oro(ig)     
     902     ENDDO
     903
     904     l2limit(1:ngrid)=0.
     905
     906     DO k=2,klev
     907        DO ig=1,ngrid
     908           winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2)
     909           IF (zlev(ig,k) .LE. h_oro(ig)) THEN  ! sous l'orographie
     910              l2strat= kapb*pbl_lmixmin_alpha*winds(ig,k)/sqrt(max(n2(ig,k),1.E-10))  ! si stratifi??, amplitude d'oscillation * kappab (voir Van de Wiel et al 2008)
     911              l2neutre=kap*zlev(ig,k)*h_oro(ig)/(kap*zlev(ig,k)+h_oro(ig))            ! Dans le cas neutre, formule de blackadar. tend asymptotiquement vers h
     912              l2neutre=MIN(l2neutre,lmax)                                             ! On majore par lmax
     913              l2limit(ig)=MIN(l2neutre,l2strat)                                       ! Calcule de l2 (minimum de la longueur en cas neutre et celle en situation stratifi??e)
     914              l2(ig,k)=l2limit(ig)
     915                                     
     916           ELSE IF (zlev(ig,k) .LE. hlim(ig)) THEN ! Si on est au dessus des montagnes, mais affect?? encore par elles
     917
     918      ! Au dessus des montagnes, on prend la l2limit au sommet des montagnes
     919      ! (la derni??re calcul??e dans la boucle k, vu que k est un indice croissant avec z)
     920      ! et on multiplie l2limit par une fonction qui d??croit entre h et hlim
     921              l2(ig,k)=l2limit(ig)*famorti(zlev(ig,k),h_oro(ig), hlim(ig))
     922           ELSE                                                                    ! Au dessus de extent*h, on prend l2=l0
     923              l2(ig,k)=0.
     924           END IF
     925        ENDDO
     926     ENDDO
     927   ENDIF                                                                        ! pbl_lmixmin_alpha
     928
     929!==================================================================================
     930! On prend le max entre la longueur de melange de blackadar et celle calcul??e
     931! en fonction de la topographie
     932!===================================================================================
     933
     934
     935 DO k=2,klev
     936    DO ig=1,ngrid
     937       lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin)
     938   ENDDO
     939 ENDDO
     940
     941! Diagnostics
     942
     943 DO k=2,klev
     944    DO ig=1,ngrid
     945       jg=ni(ig)
     946       l_mix(jg,k,nsrf)=lmix(ig,k)
     947       l_mixmin(jg,k,nsrf)=l2(ig,k)
     948    ENDDO
     949 ENDDO
     950 DO ig=1,ngrid
     951    jg=ni(ig)
     952    l_mix(jg,1,nsrf)=hlim(ig)
     953 ENDDO
     954
     955
     956
     957END SUBROUTINE mixinglength
Note: See TracChangeset for help on using the changeset viewer.