Ignore:
Timestamp:
Dec 10, 2009, 10:02:56 AM (14 years ago)
Author:
Laurent Fairhead
Message:

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

Location:
LMDZ4/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk

  • LMDZ4/trunk/libf/phylmd/physiq.F

    r1149 r1279  
    1 c
     1! $Id$
     2!
    23c#define IO_DEBUG
    34
    45      SUBROUTINE physiq (nlon,nlev,
    5      .            debut,lafin,rjourvrai,gmtime,pdtphys,
     6     .            debut,lafin,jD_cur, jH_cur,pdtphys,
    67     .            paprs,pplay,pphi,pphis,presnivs,clesphy0,
    78     .            u,v,t,qx,
    … …  
    1112     .            , PVteta)
    1213
    13       USE ioipsl
     14      USE ioipsl, only: histbeg, histvert, histdef, histend, histsync,
     15     $     histwrite, ju2ymds, ymds2ju, ioget_year_len
    1416      USE comgeomphy
     17      USE phys_cal_mod
    1518      USE write_field_phy
    1619      USE dimphy
    … …  
    2831      USE fonte_neige_mod, ONLY  : fonte_neige_get_vars
    2932      USE phys_output_mod
     33      use open_climoz_m, only: open_climoz ! ozone climatology from a file
     34      use regr_pr_av_m, only: regr_pr_av
     35      use netcdf95, only: nf95_close
     36      use mod_phys_lmdz_mpi_data, only: is_mpi_root
     37      USE aero_mod
     38      use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
     39      use conf_phys_m, only: conf_phys
     40      use radlwsw_m, only: radlwsw
    3041
    3142      IMPLICIT none
    … …  
    5061c
    5162c nlon----input-I-nombre de points horizontaux
    52 c nlev----input-I-nombre de couches verticales
     63c nlev----input-I-nombre de couches verticales, doit etre egale a klev
    5364c debut---input-L-variable logique indiquant le premier passage
    5465c lafin---input-L-variable logique indiquant le dernier passage
    55 c rjour---input-R-numero du jour de l'experience
    56 c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
     66c jD_cur       -R-jour courant a l'appel de la physique (jour julien)
     67c jH_cur       -R-heure courante a l'appel de la physique (jour julien)
    5768c pdtphys-input-R-pas d'integration pour la physique (seconde)
    5869c paprs---input-R-pression pour chaque inter-couche (en Pa)
    … …  
    104115      PARAMETER (ok_stratus=.FALSE.)
    105116c======================================================================
    106       LOGICAL, SAVE :: rnpb=.TRUE.
    107 c$OMP THREADPRIVATE(rnpb)
    108117      REAL amn, amx
    109118      INTEGER igout
    … …  
    181190      INTEGER nlon
    182191      INTEGER nlev
    183       REAL rjourvrai
    184       REAL gmtime
     192      REAL, intent(in):: jD_cur, jH_cur
     193
    185194      REAL pdtphys
    186195      LOGICAL debut, lafin
    … …  
    279288      real T2STD(klon,nlevSTD)
    280289c
    281 #include "radepsi.h"
    282290#include "radopt.h"
    283291c
    … …  
    523531      REAL clesphy0( longcles      )
    524532c
    525 c Variables quasi-arguments
    526 c
    527       REAL xjour
    528       SAVE xjour
    529 c$OMP THREADPRIVATE(xjour)
    530 c
    531 c
    532533c Variables propres a la physique
    533 c
    534 c      INTEGER radpas
    535 c      SAVE radpas                 ! frequence d'appel rayonnement
    536 ccccccccc$OMP THREADPRIVATE(radpas)
    537 c
    538 cc      INTEGER iflag_con
    539 c
    540534      INTEGER itap
    541535      SAVE itap                   ! compteur pour la physique
    … …  
    705699c Conditions aux limites
    706700c
    707       INTEGER julien
     701!
     702      REAL :: day_since_equinox
     703! Date de l'equinoxe de printemps
     704      INTEGER, parameter :: mth_eq=3, day_eq=21
     705      REAL :: jD_eq
     706
     707      LOGICAL, parameter :: new_orbit = .true.
     708
    708709c
    709710      INTEGER lmt_pas
    710711      SAVE lmt_pas                ! frequence de mise a jour
    711712c$OMP THREADPRIVATE(lmt_pas)
     713      real zmasse(klon, llm)
     714C     (column-density of mass of air in a cell, in kg m-2)
     715      real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
    712716
    713717cIM sorties
    … …  
    731735      EXTERNAL hgardfou  ! verifier les temperatures
    732736      EXTERNAL nuage     ! calculer les proprietes radiatives
    733       EXTERNAL o3cm      ! initialiser l'ozone
     737CC      EXTERNAL o3cm      ! initialiser l'ozone
    734738      EXTERNAL orbite    ! calculer l'orbite terrestre
    735       EXTERNAL ozonecm   ! prescrire l'ozone
    736739      EXTERNAL phyetat0  ! lire l'etat initial de la physique
    737740      EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
    738       EXTERNAL radlwsw   ! rayonnements solaire et infrarouge
    739741      EXTERNAL suphel    ! initialiser certaines constantes
    740742      EXTERNAL transp    ! transport total de l'eau et de l'energie
    … …  
    877879c
    878880      REAL ratqss(klon,klev),ratqsc(klon,klev)
    879       real ratqsbas,ratqshaut
    880       save ratqsbas,ratqshaut
    881 c$OMP THREADPRIVATE(ratqsbas,ratqshaut)
     881      real ratqsbas,ratqshaut,tau_ratqs
     882      save ratqsbas,ratqshaut,tau_ratqs
     883c$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
    882884      real zpt_conv(klon,klev)
    883885
    … …  
    887889      logical ok_newmicro
    888890      save ok_newmicro
     891      real ref_liq(klon,klev), ref_ice(klon,klev)
    889892c$OMP THREADPRIVATE(ok_newmicro)
    890893      save fact_cldcon,facttemps
    … …  
    972975      REAL zx_tmp_fiNC(klon,nlevSTD)
    973976c#endif
    974       REAL*8 zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
     977      REAL(KIND=8) zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
    975978      REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
    976979      REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
    … …  
    10561059      CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
    10571060      CHARACTER*40 tinst, tave, typeval
    1058 cjq   Aerosol effects (Johannes Quaas, 27/11/2003)
    1059       REAL sulfate(klon, klev) ! SO4 aerosol concentration [ug/m3]
    1060 
    10611061      REAL cldtaupi(klon,klev)  ! Cloud optical thickness for pre-industrial (pi) aerosols
    10621062
    … …  
    10671067
    10681068      ! Aerosol optical properties
    1069 
    1070       ! Aerosol optical properties by INCA model
    1071       CHARACTER*4              ::    rfname(9)
    1072       REAL aerindex(klon)       ! POLDER aerosol index
    1073      
     1069      CHARACTER*4, DIMENSION(naero_grp) :: rfname
     1070      REAL, DIMENSION(klon)          :: aerindex     ! POLDER aerosol index
     1071      REAL, DIMENSION(klon,klev)     :: mass_solu_aero    ! total mass concentration for all soluble aerosols[ug/m3]
     1072      REAL, DIMENSION(klon,klev)     :: mass_solu_aero_pi ! - " - (pre-industrial value)
     1073      INTEGER :: naero ! aerosol species
     1074
    10741075      ! Parameters
    10751076      LOGICAL ok_ade, ok_aie    ! Apply aerosol (in)direct effects or not
    … …  
    10801081                                      ! false : lecture des aerosol dans un fichier
    10811082c$OMP THREADPRIVATE(aerosol_couple)   
    1082 
     1083      INTEGER, SAVE :: flag_aerosol
     1084c$OMP THREADPRIVATE(flag_aerosol)
     1085      LOGICAL, SAVE :: new_aod
     1086c$OMP THREADPRIVATE(new_aod)
     1087   
    10831088c
    10841089c Declaration des constantes et des fonctions thermodynamiques
    … …  
    10861091      LOGICAL,SAVE :: first=.true.
    10871092c$OMP THREADPRIVATE(first)
     1093
     1094      integer iunit
     1095
     1096      integer, save::  read_climoz ! read ozone climatology
     1097C     Allowed values are 0, 1 and 2
     1098C     0: do not read an ozone climatology
     1099C     1: read a single ozone climatology that will be used day and night
     1100C     2: read two ozone climatologies, the average day and night
     1101C     climatology and the daylight climatology
     1102
     1103      integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
     1104
     1105      real, pointer, save:: press_climoz(:)
     1106!     edges of pressure intervals for ozone climatologies, in Pa, in strictly
     1107!     ascending order
     1108
     1109      integer, save:: co3i = 0
     1110!     time index in NetCDF file of current ozone fields
     1111c$OMP THREADPRIVATE(co3i)
     1112
     1113      integer ro3i
     1114!     required time index in NetCDF file for the ozone fields, between 1
     1115!     and 360
     1116
    10881117#include "YOMCST.h"
    10891118#include "YOETHF.h"
    … …  
    10961125cIM 100106 END : pouvoir sortir les ctes de la physique
    10971126c
     1127!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1128c Declarations pour Simulateur COSP
     1129c============================================================
     1130      real :: mr_ozone(klon,klev)
    10981131c======================================================================
    10991132! Ecriture eventuelle d'un profil verticale en entree de la physique.
    … …  
    11061139         write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
    11071140         write(lunout,*)
    1108      s 'nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys'
     1141     s 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
    11091142         write(lunout,*)
    1110      s  nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys
    1111 
    1112          write(lunout,*) 'papers, play, phi, u, v, t, omega'
    1113          do k=1,nlev
     1143     s  nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
     1144
     1145         write(lunout,*) 'paprs, play, phi, u, v, t'
     1146         do k=1,klev
    11141147            write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),
    1115      s   u(igout,k),v(igout,k),t(igout,k),omega(igout,k)
     1148     s   u(igout,k),v(igout,k),t(igout,k)
    11161149         enddo
    11171150         write(lunout,*) 'ovap (g/kg),  oliq (g/kg)'
    1118          do k=1,nlev
     1151         do k=1,klev
    11191152            write(lunout,*) qx(igout,k,1)*1000,qx(igout,k,2)*1000.
    11201153         enddo
    … …  
    11321165
    11331166      torsfc=0.
     1167      forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
    11341168
    11351169      if (first) then
    … …  
    11401174      print*, 'Allocation des variables locales et sauvegardees'
    11411175      call phys_local_var_init
    1142       call phys_state_var_init
     1176c     appel a la lecture du run.def physique
     1177      call conf_phys(ok_journe, ok_mensuel,
     1178     .     ok_instan, ok_hf,
     1179     .     ok_LES,
     1180     .     solarlong0,seuil_inversion,
     1181     .     fact_cldcon, facttemps,ok_newmicro,iflag_radia,
     1182     .     iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,
     1183     .     ok_ade, ok_aie, aerosol_couple,
     1184     .     flag_aerosol, new_aod,
     1185     .     bl95_b0, bl95_b1,
     1186     .     iflag_thermals,nsplit_thermals,tau_thermals,
     1187     .     iflag_thermals_ed,iflag_thermals_optflux,
     1188c     nv flags pour la convection et les poches froides
     1189     .     iflag_coupl,iflag_clos,iflag_wake, read_climoz)
     1190      call phys_state_var_init(read_climoz)
    11431191      print*, '================================================='
    11441192
    … …  
    11561204        first=.false.
    11571205
    1158       endif  ! fisrt
     1206      endif  ! first
    11591207
    11601208       modname = 'physiq'
    … …  
    11751223
    11761224c======================================================================
    1177       xjour = rjourvrai
     1225! Gestion calendrier : mise a jour du module phys_cal_mod
     1226!
     1227      CALL phys_cal_update(jD_cur,jH_cur)
     1228
    11781229c
    11791230c Si c'est le debut, il faut initialiser plusieurs choses
    … …  
    11811232c
    11821233       IF (debut) THEN
    1183 C
    11841234!rv
    11851235cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
    … …  
    11901240         u10m(:,:)=0.
    11911241         v10m(:,:)=0.
    1192          piz_ae(:,:,:)=0.
    1193          tau_ae(:,:,:)=0.
    1194          cg_ae(:,:,:)=0.
    11951242         rain_con(:)=0.
    11961243         snow_con(:)=0.
    1197          bl95_b0=0.
    1198          bl95_b1=0.
    11991244         topswai(:)=0.
    12001245         topswad(:)=0.
    … …  
    12051250         wmax_th(:)=0.
    12061251         tau_overturning_th(:)=0.
     1252
    12071253         IF (config_inca /= 'none') THEN
    1208             tau_inca(:,:,:,:) = 0.
    1209             piz_inca(:,:,:,:) = 0.
    1210             cg_inca(:,:,:,:)  = 0.
    1211             ccm(:,:,:)        = 0.
    1212             topswai_inca(:)   = 0.
    1213             topswad_inca(:)   = 0.
    1214             topswad0_inca(:)  = 0.
    1215             topsw_inca(:,:)   = 0.
    1216             topsw0_inca(:,:)  = 0.
    1217             solswai_inca(:)   = 0.
    1218             solswad_inca(:)   = 0.
    1219             solswad0_inca(:)  = 0.
    1220             solsw_inca(:,:)   = 0.
    1221             solsw0_inca(:,:)  = 0.
     1254            ! jg : initialisation jusqu'au ces variables sont dans restart
     1255            ccm(:,:,:) = 0.
     1256            tau_aero(:,:,:,:) = 0.
     1257            piz_aero(:,:,:,:) = 0.
     1258            cg_aero(:,:,:,:) = 0.
    12221259         END IF
    12231260
    … …  
    12301267         IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
    12311268c
    1232 c appel a la lecture du run.def physique
    1233 c
    1234          call conf_phys(ok_journe, ok_mensuel,
    1235      .                  ok_instan, ok_hf,
    1236      .                  ok_LES,
    1237      .                  solarlong0,seuil_inversion,
    1238      .                  fact_cldcon, facttemps,ok_newmicro,iflag_radia,
    1239      .                  iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,
    1240      .                  ok_ade, ok_aie, aerosol_couple,
    1241      .                  bl95_b0, bl95_b1,
    1242      .                  iflag_thermals,nsplit_thermals,tau_thermals,
    1243      .                  iflag_thermals_ed,iflag_thermals_optflux,
    1244 cnv flags pour la convection et les poches froides
    1245      .                   iflag_coupl,iflag_clos,iflag_wake)
    1246 
    12471269      print*,'iflag_coupl,iflag_clos,iflag_wake',
    12481270     .   iflag_coupl,iflag_clos,iflag_wake
    … …  
    13771399         ENDIF
    13781400
    1379            rugoro=0.
     1401           DO i=1,klon
     1402             rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
     1403           ENDDO
     1404
    13801405c34EK
    13811406         IF (ok_orodr) THEN
    1382 
    1383            rugoro=0.
    13841407
    13851408!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    14131436     .                   lmt_pas
    14141437c
    1415 cIM200505        ecrit_mth = NINT(86400./dtime *ecritphy)  ! tous les ecritphy jours
    1416 c        IF (ok_mensuel) THEN
    1417 c        WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
    1418 c    .                   ecrit_mth
    1419 c        ENDIF
    1420 c        ecrit_day = NINT(86400./dtime *1.0)  ! tous les jours
    1421 c        IF (ok_journe) THEN
    1422 c        WRITE(lunout,*)'La frequence de sortie journaliere est de ',
    1423 c    .                   ecrit_day
    1424 c        ENDIF
    1425 cIM 130904 BEG
    1426 cIM 080205      ecrit_hf = 86400./dtime *0.25  ! toutes les 6h
    1427 cIM 170305     
    1428 c        ecrit_hf = 86400./dtime/12.  ! toutes les 2h
    1429 cIM 230305     
    1430 cIM200505        ecrit_hf = 86400./dtime *0.25  ! toutes les 6h
    1431 c
    1432 cIM200505        ecrit_hf2mth = ecrit_day/ecrit_hf*30
    1433 c
    1434 cIM200505        IF (ok_journe) THEN
    1435 cIM200505        WRITE(lunout,*)'La frequence de sortie hf est de ',
    1436 cIM200505    .                   ecrit_hf
    1437 cIM200505        ENDIF
    1438 cIM 130904 END
    1439 ccc         ecrit_ins = NINT(86400./dtime *0.5)  ! 2 fois par jour
    1440 ccc         ecrit_ins = NINT(86400./dtime *0.25)  ! 4 fois par jour
    1441 c        ecrit_ins = NINT(86400./dtime/48.)  ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
    1442 c        ecrit_ins = NINT(86400./dtime/12.)  ! toutes les deux heures
    1443 cIM200505        ecrit_ins = NINT(86400./dtime/8.)  ! toutes les trois heures
    1444 cIM200505        IF (ok_instan) THEN
    1445 cIM200505        WRITE(lunout,*)'La frequence de sortie instant. est de ',
    1446 cIM200505    .                   ecrit_ins
    1447 cIM200505        ENDIF
    1448 cIM200505        ecrit_reg = NINT(86400./dtime *0.25)  ! 4 fois par jour
    1449 cIM200505        IF (ok_region) THEN
    1450 cIM200505        WRITE(lunout,*)'La frequence de sortie region est de ',
    1451 cIM200505    .                   ecrit_reg
    1452 cIM200505        ENDIF
    1453 cIM 030306 BEG
    1454 cIM ecrit_hf2mth = nombre de pas de temps de calcul de hf par mois apres lequel on ecrit
    1455 cIM : ne pas modifier ecrit_hf2mth
    1456 c
    1457 cIM 250308bad guide        ecrit_hf2mth = 30*1/ecrit_hf
    1458          ecrit_hf2mth = ecrit_mth/ecrit_hf
    1459 c ecrit_ins en secondes, chaque pas de temps de la physique
    1460          ecrit_ins = dtime
    1461 cIM on passe les frequences de jours en secondes : ecrit_ins, ecrit_hf, ecrit_day, ecrit_mth, ecrit_tra, ecrit_reg
    1462          ecrit_hf = ecrit_hf * un_jour
    1463 !IM
    1464          IF(ecrit_day.LE.1.) THEN
    1465           ecrit_day = ecrit_day * un_jour !en secondes
    1466          ENDIF
    1467 !IM
    1468          ecrit_mth = ecrit_mth * un_jour
    1469          ecrit_reg = ecrit_reg * un_jour
    1470          ecrit_tra = ecrit_tra * un_jour
    1471          ecrit_ISCCP = ecrit_ISCCP * un_jour
    1472          ecrit_LES = ecrit_LES * un_jour
    1473 c
    1474          PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
    1475      .   ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
    1476      .   ecrit_hf2mth
    14771438cIM 030306 END
    14781439
    … …  
    14941455c$OMP MASTER
    14951456       call phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta,
    1496      &                        ctetaSTD,dtime,presnivs,ok_veget,
     1457     &                        ctetaSTD,dtime,ok_veget,
    14971458     &                        type_ocean,iflag_pbl,ok_mensuel,ok_journe,
    1498      &                        ok_hf,ok_instan,ok_LES,ok_ade,ok_aie)
     1459     &                        ok_hf,ok_instan,ok_LES,ok_ade,ok_aie,
     1460     &                        read_climoz, new_aod, aerosol_couple)
    14991461c$OMP END MASTER
    15001462c$OMP BARRIER
    … …  
    15141476#endif
    15151477
     1478cIM 250308bad guide        ecrit_hf2mth = 30*1/ecrit_hf
     1479         ecrit_hf2mth = ecrit_mth/ecrit_hf
     1480
     1481         ecrit_hf = ecrit_hf * un_jour
     1482!IM
     1483         IF(ecrit_day.LE.1.) THEN
     1484          ecrit_day = ecrit_day * un_jour !en secondes
     1485         ENDIF
     1486!IM
     1487         ecrit_mth = ecrit_mth * un_jour
     1488         ecrit_ins = ecrit_ins * un_jour
     1489         ecrit_reg = ecrit_reg * un_jour
     1490         ecrit_tra = ecrit_tra * un_jour
     1491         ecrit_ISCCP = ecrit_ISCCP * un_jour
     1492         ecrit_LES = ecrit_LES * un_jour
     1493c
     1494         PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
     1495     .   ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
     1496     .   ecrit_hf2mth
     1497cIM 030306 END
     1498
     1499
    15161500cXXXPB Positionner date0 pour initialisation de ORCHIDEE
    1517       date0 = zjulian
    1518 C      date0 = day_ini
     1501      date0 = jD_ref
    15191502      WRITE(*,*) 'physiq date0 : ',date0
    15201503c
    … …  
    15341517         CALL VTe(VTphysiq)
    15351518         CALL VTb(VTinca)
    1536          iii = MOD(NINT(xjour),360)
    1537          calday = FLOAT(iii) + gmtime
    1538          WRITE(lunout,*) 'initial time ', xjour, calday
     1519!         iii = MOD(NINT(xjour),360)
     1520!         calday = FLOAT(iii) + jH_cur
     1521         calday = FLOAT(days_elapsed) + jH_cur
     1522         WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
    15391523
    15401524         CALL chemini(
    … …  
    15501534     $                   pdtphys,
    15511535     $                   annee_ref,
    1552      $                   day_ini)
     1536     $                   day_ref,
     1537     $                   itau_phy)
    15531538
    15541539         CALL VTe(VTinca)
    … …  
    15641549      call iniradia(klon,klev,paprs(1,1:klev+1))
    15651550
     1551C$omp single
     1552      if (read_climoz >= 1) then
     1553         call open_climoz(ncid_climoz, press_climoz)
     1554      END IF
     1555C$omp end single
    15661556      ENDIF
    15671557!
    … …  
    15721562!
    15731563      itap   = itap + 1
    1574       julien = MOD(NINT(xjour),360)
    1575       if (julien .eq. 0) julien = 360
    1576 
    15771564!
    15781565! Update fraction of the sub-surfaces (pctsrf) and
    … …  
    15801567! on the surface fraction.
    15811568!
    1582       CALL change_srf_frac(itap, dtime, julien,
     1569      CALL change_srf_frac(itap, dtime, days_elapsed+1,
    15831570     *     pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
    1584 
    15851571
    15861572! Tendances bidons pour les processus qui n'affectent pas certaines
    … …  
    17311717c Prescrire l'ozone et calculer l'albedo sur l'ocean.
    17321718c
    1733       IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
    1734          if(prt_level.ge.1) WRITE(lunout,*)' PHYS cond  julien ',julien
    1735          CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
     1719      if (read_climoz >= 1) then
     1720C        Ozone from a file
     1721!        Update required ozone index:
     1722         ro3i = int((days_elapsed + jh_cur - jh_1jan)
     1723     $        / ioget_year_len(year_cur) * 360.) + 1
     1724         if (ro3i == 361) ro3i = 360
     1725C        (This should never occur, except perhaps because of roundup
     1726C        error. See documentation.)
     1727         if (ro3i /= co3i) then
     1728C           Update ozone field:
     1729            if (read_climoz == 1) then
     1730               call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i,
     1731     $              press_in_edg=press_climoz, paprs=paprs, v3=wo)
     1732            else
     1733C              read_climoz == 2
     1734               call regr_pr_av(ncid_climoz,
     1735     $              (/"tro3         ", "tro3_daylight"/),
     1736     $              julien=ro3i, press_in_edg=press_climoz, paprs=paprs,
     1737     $              v3=wo)
     1738            end if
     1739!           Convert from mole fraction of ozone to column density of ozone in a
     1740!           cell, in kDU:
     1741            forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l)
     1742     $           * rmo3 / rmd * zmasse / dobson_u / 1e3
     1743C           (By regridding ozone values for LMDZ only once every 360th of
     1744C           year, we have already neglected the variation of pressure in one
     1745C           360th of year. So do not recompute "wo" at each time step even if
     1746C           "zmasse" changes a little.)
     1747            co3i = ro3i
     1748         end if
     1749      elseif (MOD(itap-1,lmt_pas) == 0) THEN
     1750C        Once per day, update ozone from Royer:
     1751         wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
    17361752      ENDIF
    17371753c
    … …  
    17741790! doit donc etre placé avant radlwsw et pbl_surface
    17751791
     1792! calcul selon la routine utilisee pour les planetes
     1793      if (new_orbit) then
     1794        call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
     1795        day_since_equinox = (jD_cur + jH_cur) - jD_eq
     1796!        day_since_equinox = (jD_cur) - jD_eq
     1797        call solarlong(day_since_equinox, zlongi, dist)
     1798      else     
     1799! calcul selon la routine utilisee pour l'AR4
    17761800!   choix entre calcul de la longitude solaire vraie ou valeur fixee a
    17771801!   solarlong0
    1778 
    1779       if (solarlong0<-999.) then
    1780          CALL orbite(FLOAT(julien),zlongi,dist)
    1781       else
    1782          zlongi=solarlong0  ! longitude solaire vraie
    1783          dist=1.            ! distance au soleil / moyenne
     1802        if (solarlong0<-999.) then
     1803           CALL orbite(FLOAT(days_elapsed+1),zlongi,dist)
     1804        else
     1805           zlongi=solarlong0  ! longitude solaire vraie
     1806           dist=1.            ! distance au soleil / moyenne
     1807        endif
    17841808      endif
    1785 
    1786       if(prt_level.ge.1) print*,'Longitude solaire ',zlongi,solarlong0
     1809      if(prt_level.ge.1)                                                &
     1810     &    write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
    17871811
    17881812!  Avec ou sans cycle diurne
    17891813      IF (cycle_diurne) THEN
    17901814        zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
    1791         CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
     1815        CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
    17921816      ELSE
    17931817        CALL angle(zlongi, rlat, fract, rmu0)
    … …  
    18221846
    18231847      CALL pbl_surface(
    1824      e     dtime,     date0,     itap,    julien,
     1848     e     dtime,     date0,     itap,    days_elapsed+1,
    18251849     e     debut,     lafin,
    18261850     e     rlon,      rlat,      rugoro,  rmu0,     
    … …  
    19331957         END DO
    19341958      END DO
     1959      if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ',
     1960     $     omega(igout, :)
    19351961
    19361962      IF (iflag_con.EQ.1) THEN
    … …  
    21192145
    21202146            if (itop_con(i).gt.klev-3) then
    2121                print*,'La convection monte trop haut '
    2122                print*,'itop_con(,',i,',)=',itop_con(i)
     2147              if(prt_level >= 9) then
     2148                write(lunout,*)'La convection monte trop haut '
     2149                write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
     2150              endif
    21232151            endif
    21242152          ENDDO
    … …  
    25192547         enddo
    25202548
    2521 !   les ratqs sont une conbinaison de ratqss et ratqsc
    2522 !   1800s, en dur pour le moment, est le temps de
    2523 !   relaxation des ratqs
    2524 
    2525          facteur=exp(-pdtphys/1800.)
    2526 
    2527          print*,'WARNING ratqs a revoir '
     2549!   les ratqs sont une combinaison de ratqss et ratqsc
     2550       print*,'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
     2551
     2552         if (tau_ratqs>1.e-10) then
     2553            facteur=exp(-pdtphys/tau_ratqs)
     2554         else
     2555            facteur=0.
     2556         endif
    25282557         ratqs(:,:)=ratqsc(:,:)*(1.-facteur)+ratqs(:,:)*facteur
    2529          ratqs(:,:)=sqrt(ratqs(:,:)**2+ratqss(:,:)**2)
    2530 
     2558!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2559! FH 22/09/2009
     2560! La ligne ci-dessous faisait osciller le modele et donnait une solution
     2561! assymptotique bidon et dépendant fortement du pas de temps.
     2562!        ratqs(:,:)=sqrt(ratqs(:,:)**2+ratqss(:,:)**2)
     2563!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2564         ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
    25312565      else
    25322566!   on ne prend que le ratqs stable pour fisrtilp
    … …  
    26582692cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
    26592693      IF (ok_ade.OR.ok_aie) THEN
    2660        IF ( .NOT. aerosol_couple ) THEN
    2661          ! Get sulfate aerosol distribution
    2662          CALL readsulfate(rjourvrai, debut, sulfate)
    2663          CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
    2664 
    2665          ! Calculate aerosol optical properties (Olivier Boucher)
    2666          CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
    2667      .        tau_ae, piz_ae, cg_ae, aerindex)
    2668        ENDIF
     2694         IF (.NOT. aerosol_couple)
     2695     &        CALL readaerosol_optic(
     2696     &        debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref,
     2697     &        pdtphys, pplay, paprs, t_seri, rhcl, presnivs,
     2698     &        mass_solu_aero, mass_solu_aero_pi,
     2699     &        tau_aero, piz_aero, cg_aero,
     2700     &        tausum_aero, tau3d_aero)
    26692701      ELSE
    2670         tau_ae(:,:,:)=0.0
    2671         piz_ae(:,:,:)=0.0
    2672         cg_ae(:,:,:)=0.0
     2702         tau_aero(:,:,:,:) = 0.
     2703         piz_aero(:,:,:,:) = 0.
     2704         cg_aero(:,:,:,:)  = 0.
    26732705      ENDIF
    26742706
    … …  
    27912823         CALL VTe(VTphysiq)
    27922824         CALL VTb(VTinca)
    2793          calday = FLOAT(julien) + gmtime
     2825         calday = FLOAT(days_elapsed + 1) + jH_cur
    27942826
    27952827         IF (config_inca == 'aero') THEN
    27962828            CALL AEROSOL_METEO_CALC(
    27972829     $           calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
    2798      $           prfl,psfl,pctsrf,airephy,xjour,rlat,rlon,u10m,v10m)
     2830     $           prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
    27992831         END IF
    28002832
    … …  
    28022834
    28032835         CALL chemhook_begin (calday,
    2804      $                          julien,
    2805      $                          gmtime,
     2836     $                          days_elapsed+1,
     2837     $                          jH_cur,
    28062838     $                          pctsrf(1,1),
    28072839     $                          rlat,
    … …  
    28152847     $                          u,
    28162848     $                          v,
    2817      $                          wo,
     2849     $                          wo(:, :, 1),
    28182850     $                          q_seri,
    28192851     $                          zxtsol,
    … …  
    28472879
    28482880      IF (aerosol_couple) THEN
    2849          sulfate(:,:) = ccm(:,:,1)
    2850          sulfate_pi(:,:) = ccm(:,:,2)
    2851       ENDIF
     2881         mass_solu_aero(:,:)    = ccm(:,:,1)
     2882         mass_solu_aero_pi(:,:) = ccm(:,:,2)
     2883      END IF
    28522884
    28532885      if (ok_newmicro) then
    … …  
    28572889     .            flwp, fiwp, flwc, fiwc,
    28582890     e            ok_aie,
    2859      e            sulfate, sulfate_pi,
     2891     e            mass_solu_aero, mass_solu_aero_pi,
    28602892     e            bl95_b0, bl95_b1,
    2861      s            cldtaupi, re, fl)
     2893     s            cldtaupi, re, fl, ref_liq, ref_ice)
    28622894      else
    28632895      CALL nuage (paprs, pplay,
    … …  
    28652897     .            cldh, cldl, cldm, cldt, cldq,
    28662898     e            ok_aie,
    2867      e            sulfate, sulfate_pi,
     2899     e            mass_solu_aero, mass_solu_aero_pi,
    28682900     e            bl95_b0, bl95_b1,
    28692901     s            cldtaupi, re, fl)
    … …  
    28952927      IF (aerosol_couple) THEN
    28962928#ifdef INCA
    2897       CALL radlwsw_inca
    2898      e            (kdlon,kflev,dist, rmu0, fract, solaire,
    2899      e             paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
    2900      e             wo,
    2901      e             cldfra, cldemi, cldtau,
    2902      s             heat,heat0,cool,cool0,radsol,albpla,
    2903      s             topsw,toplw,solsw,sollw,
    2904      s             sollwdown,
    2905      s             topsw0,toplw0,solsw0,sollw0,
    2906      s             lwdn0, lwdn, lwup0, lwup,
    2907      s             swdn0, swdn, swup0, swup,
    2908      e             ok_ade, ok_aie,
    2909      e             tau_inca, piz_inca, cg_inca,
    2910      s             topswad_inca, solswad_inca,
    2911      s             topswad0_inca, solswad0_inca,
    2912      s             topsw_inca, topsw0_inca,
    2913      s             solsw_inca, solsw0_inca,
    2914      e             cldtaupi,
    2915      s             topswai_inca, solswai_inca)
     2929         CALL radlwsw_inca
     2930     e        (kdlon,kflev,dist, rmu0, fract, solaire,
     2931     e        paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
     2932     e        wo(:, :, 1),
     2933     e        cldfra, cldemi, cldtau,
     2934     s        heat,heat0,cool,cool0,radsol,albpla,
     2935     s        topsw,toplw,solsw,sollw,
     2936     s        sollwdown,
     2937     s        topsw0,toplw0,solsw0,sollw0,
     2938     s        lwdn0, lwdn, lwup0, lwup,
     2939     s        swdn0, swdn, swup0, swup,
     2940     e        ok_ade, ok_aie,
     2941     e        tau_aero, piz_aero, cg_aero,
     2942     s        topswad_aero, solswad_aero,
     2943     s        topswad0_aero, solswad0_aero,
     2944     s        topsw_aero, topsw0_aero,
     2945     s        solsw_aero, solsw0_aero,
     2946     e        cldtaupi,
     2947     s        topswai_aero, solswai_aero)
     2948           
    29162949#endif
    29172950      ELSE
    2918       CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
    2919      e            (dist, rmu0, fract,
    2920      e             paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
    2921      e             wo,
    2922      e             cldfra, cldemi, cldtau,
    2923      s             heat,heat0,cool,cool0,radsol,albpla,
    2924      s             topsw,toplw,solsw,sollw,
    2925      s             sollwdown,
    2926      s             topsw0,toplw0,solsw0,sollw0,
    2927      s             lwdn0, lwdn, lwup0, lwup,
    2928      s             swdn0, swdn, swup0, swup,
    2929      e             ok_ade, ok_aie, ! new for aerosol radiative effects
    2930      e             tau_ae, piz_ae, cg_ae, ! ="=
    2931      s             topswad, solswad, ! ="=
    2932      e             cldtaupi, ! ="=
    2933      s             topswai, solswai,zqsat,flwc,fiwc) ! ="=
    2934       ENDIF
     2951
     2952         CALL radlwsw
     2953     e        (dist, rmu0, fract,
     2954     e        paprs, pplay,zxtsol,albsol1, albsol2,
     2955     e        t_seri,q_seri,wo,
     2956     e        cldfra, cldemi, cldtau,
     2957     e        ok_ade, ok_aie,
     2958     e        tau_aero, piz_aero, cg_aero,
     2959     e        cldtaupi,new_aod,
     2960     e        zqsat, flwc, fiwc,
     2961     s        heat,heat0,cool,cool0,radsol,albpla,
     2962     s        topsw,toplw,solsw,sollw,
     2963     s        sollwdown,
     2964     s        topsw0,toplw0,solsw0,sollw0,
     2965     s        lwdn0, lwdn, lwup0, lwup,
     2966     s        swdn0, swdn, swup0, swup,
     2967     s        topswad_aero, solswad_aero,
     2968     s        topswai_aero, solswai_aero,
     2969     o        topswad0_aero, solswad0_aero,
     2970     o        topsw_aero, topsw0_aero,
     2971     o        solsw_aero, solsw0_aero,
     2972     o        topswcf_aero, solswcf_aero)
     2973         
     2974
     2975      ENDIF ! aerosol_couple
    29352976      itaprad = 0
    2936       ENDIF
     2977      ENDIF ! MOD(itaprad,radpas)
    29372978      itaprad = itaprad + 1
    29382979
    … …  
    31233164cIM calcul composantes axiales du moment angulaire et couple des montagnes
    31243165c
    3125       IF (is_sequential .AND. ok_orodr .AND. ok_orolf ) THEN
     3166      IF (is_sequential) THEN
    31263167     
    3127         CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
     3168        CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
    31283169     C                 ra,rg,romega,
    31293170     C                 rlat,rlon,pphis,
    … …  
    31433184c
    31443185c
     3186!====================================================================
     3187! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
     3188!====================================================================
     3189! Abderrahmane 24.08.09
     3190
     3191      IF (ok_cosp) THEN
     3192! adeclarer
     3193#ifdef CPP_COSP
     3194       IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
     3195
     3196       print*,'freq_cosp',freq_cosp
     3197          mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
     3198!       print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
     3199!     s        ref_liq,ref_ice
     3200          call phys_cosp(itap,dtime,freq_cosp,
     3201     $                 ecrit_mth,ecrit_day,ecrit_hf,overlap,
     3202     $                   klon,klev,rlon,rlat,presnivs,
     3203     $                   ref_liq,ref_ice,
     3204     $                   pctsrf(:,is_ter)+pctsrf(:,is_lic),
     3205     $                   zu10m,zv10m,
     3206     $                   zphi,paprs(:,1:klev),pplay,zxtsol,t_seri,
     3207     $                   qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc,
     3208     $                   prfl(:,1:klev),psfl(:,1:klev),
     3209     $                   pmflxr(:,1:klev),pmflxs(:,1:klev),
     3210     $                   mr_ozone,cldtau, cldemi)
     3211!     L          calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
     3212!     L          cfaddbze,clcalipso2,dbze,cltlidarradar,
     3213!     M          clMISR,
     3214!     R          clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
     3215!     I          tauisccp,albisccp,meantbisccp,meantbclrisccp)
     3216
     3217         ENDIF
     3218
     3219#endif
     3220       ENDIF  !ok_cosp
     3221!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    31453222cAA
    31463223cAA Installation de l'interface online-offline pour traceurs
    … …  
    31503227c====================================================================
    31513228C
    3152       IF (config_inca /= 'none') rnpb=.FALSE.
    3153 
    3154       call phytrac (     rnpb,
    3155      I                   itap,
    3156      I                   julien,
    3157      I                   gmtime,
    3158      I                   debut,
    3159      I                   lafin,
    3160      I                   nlon,
    3161      I                   nlev,
    3162      I                   dtime,
    3163      I                   u,
    3164      I                   v,
    3165      I                   t,
    3166      I                   paprs,
    3167      I                   pplay,
    3168      I                   pmfu,
    3169      I                   pmfd,
    3170      I                   pen_u,
    3171      I                   pde_u,
    3172      I                   pen_d,
    3173      I                   pde_d,
    3174      I                   cdragh,
    3175      I                   coefh,
    3176      I                   fm_therm,
    3177      I                   entr_therm,
    3178      I                   u1,
    3179      I                   v1,
    3180      I                   ftsol,
    3181      I                   pctsrf,
    3182      I                   rlat,
    3183      I                   frac_impa,
    3184      I                   frac_nucl,
    3185      I                   rlon,
    3186      I                   presnivs,
    3187      I                   pphis,
    3188      I                   pphi,
    3189      I                   albsol1,
    3190      I                   qx(1,1,1),
    3191      I                   rhcl,
    3192      I                   cldfra,
    3193      I                   rneb,
    3194      I                   diafra,
    3195      I                   cldliq,
    3196      I                   itop_con,
    3197      I                   ibas_con,
    3198      I                   pmflxr,
    3199      I                   pmflxs,
    3200      I                   prfl,
    3201      I                   psfl,
    3202      I                   da,
    3203      I                   phi,
    3204      I                   mp,
    3205      I                   upwd,
    3206      I                   dnwd,
    3207      I                   aerosol_couple,
    3208      I                   flxmass_w,
    3209      I                   tau_inca,
    3210      I                   piz_inca,
    3211      I                   cg_inca,
    3212      I                   ccm,
    3213      I                   rfname,
    3214      O                   tr_seri)
     3229
     3230      call phytrac (
     3231     I     itap,     days_elapsed+1,    jH_cur,   debut,
     3232     I     lafin,    dtime,     u, v,     t,
     3233     I     paprs,    pplay,     pmfu,     pmfd,
     3234     I     pen_u,    pde_u,     pen_d,    pde_d,
     3235     I     cdragh,   coefh,     fm_therm, entr_therm,
     3236     I     u1,       v1,        ftsol,    pctsrf,
     3237     I     rlat,     frac_impa, frac_nucl,rlon,
     3238     I     presnivs, pphis,     pphi,     albsol1,
     3239     I     qx(:,:,ivap),rhcl,   cldfra,   rneb,
     3240     I     diafra,   cldliq,    itop_con, ibas_con,
     3241     I     pmflxr,   pmflxs,    prfl,     psfl,
     3242     I     da,       phi,       mp,       upwd,     
     3243     I     dnwd,     aerosol_couple,      flxmass_w,
     3244     I     tau_aero, piz_aero,  cg_aero,  ccm,
     3245     I     rfname,
     3246     O     tr_seri)
    32153247
    32163248      IF (offline) THEN
    … …  
    32183250         print*,'Attention on met a 0 les thermiques pour phystoke'
    32193251         call phystokenc (
    3220      I                   nlon,nlev,pdtphys,rlon,rlat,
     3252     I                   nlon,klev,pdtphys,rlon,rlat,
    32213253     I                   t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
    32223254     I                   fm_therm,entr_therm,
    … …  
    32523284        d_t_ec(i,k)=0.5/ZRCPD
    32533285     $      *(u(i,k)**2+v(i,k)**2-u_seri(i,k)**2-v_seri(i,k)**2)
     3286      ENDDO
     3287      ENDDO
     3288
     3289      DO k = 1, klev
     3290      DO i = 1, klon
    32543291        t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
    32553292        d_t_ec(i,k) = d_t_ec(i,k)/dtime
    … …  
    32673304C     est egale a la variation de la physique au pas de temps precedent.
    32683305C     Donc la somme de ces 2 variations devrait etre nulle.
     3306
    32693307        call diagphy(airephy,ztit,ip_ebil_phy
    32703308     e      , topsw, toplw, solsw, sollw, sens
    … …  
    33493387     $                        day_ini,
    33503388     $                        airephy,
    3351      $                        xjour,
    33523389     $                        pphi,
    33533390     $                        pphis,
    … …  
    34153452      write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
    34163453      write(lunout,*)
    3417      s 'nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos'
     3454     s 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
    34183455      write(lunout,*)
    3419      s  nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys,
     3456     s  nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,
    34203457     s  pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
    34213458     s  pctsrf(igout,is_sic)
    34223459      write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
    3423       do k=1,nlev
     3460      do k=1,klev
    34243461         write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
    34253462     s   d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
    … …  
    34273464      enddo
    34283465      write(lunout,*) 'cool,heat'
    3429       do k=1,nlev
     3466      do k=1,klev
    34303467         write(lunout,*) cool(igout,k),heat(igout,k)
    34313468      enddo
    34323469
    34333470      write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
    3434       do k=1,nlev
     3471      do k=1,klev
    34353472         write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
    34363473     s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
    … …  
    34393476      write(lunout,*) 'd_ps ',d_ps(igout)
    34403477      write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
    3441       do k=1,nlev
     3478      do k=1,klev
    34423479         write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
    34433480     s  d_qx(igout,k,1),d_qx(igout,k,2)
    … …  
    35013538!         write(97) u_seri,v_seri,t_seri,q_seri
    35023539!         close(97)
     3540C$OMP MASTER
     3541         if (read_climoz >= 1) then
     3542            if (is_mpi_root) then
     3543               call nf95_close(ncid_climoz)
     3544            end if
     3545            deallocate(press_climoz) ! pointer
     3546         end if
     3547C$OMP END MASTER
    35033548      ENDIF
    35043549     
     3550!      first=.false.
    35053551
    35063552      RETURN
Note: See TracChangeset for help on using the changeset viewer.