Changeset 1130


Ignore:
Timestamp:
Dec 20, 2013, 4:04:56 PM (11 years ago)
Author:
emillour
Message:

Mars GCM:
Series of changes to enable running in parallel (using LMDZ.COMMON dynamics);
Current LMDZ.MARS can still notheless be compiled and run in serial mode
"as previously".
Summary of main changes:

  • Main programs (newstart, start2archive, xvik) that used to be in dyn3d have been moved to phymars.
  • dyn3d/control.h is now module control_mod.F90
  • rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
  • added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallel case)
  • created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the min and max of a field over the whole planet.

EM

Location:
trunk/LMDZ.MARS
Files:
17 added
5 deleted
35 edited
3 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r1127 r1130  
    19651965- Bug fix in thermcall_main_mars, layer thickness "zdz" must be recomputed
    19661966  in every do ig=1,ngrid loop.
     1967
     1968== 20/12/2013 == EM
     1969Series of changes to enable running in parallel (using LMDZ.COMMON dynamics);
     1970Current LMDZ.MARS can still notheless be compiled and run in serial mode
     1971"as previously".
     1972Summary of main changes:
     1973- Main programs (newstart, start2archive, xvik) that used to be in dyn3d have
     1974  been moved to phymars.
     1975- dyn3d/control.h is now module control_mod.F90
     1976- rearanged input/outputs routines everywhere to handle serial/MPI cases.
     1977  physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write
     1978  routines for startfi files are gathered in module iostart.F90
     1979- added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90,
     1980  dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90,
     1981  mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90,
     1982  mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90,
     1983  mod_phys_lmdz_transfert_para.F90 in phymars
     1984  and mod_const_mpi.F90 in dyn3d (for compliance with parallel case)
     1985- created generic routines 'planetwide_maxval' and 'planetwide_minval', in
     1986  module "planetwide_mod", that enable obtaining the min and max of a field
     1987  over the whole planet.

  • trunk/LMDZ.MARS/libf/dyn3d/calfis.F

    r1036 r1130  
    7171#include "comvert.h"
    7272#include "comgeom2.h"
    73 #include "control.h"
     73!#include "control.h"
    7474
    7575c    Arguments :
    … …  
    167167
    168168c
    169       IF (firstcal) THEN
    170          latfi(1)=rlatu(1)
    171          lonfi(1)=0.
    172          DO j=2,jjm
    173             DO i=1,iim
    174                latfi((j-2)*iim+1+i)= rlatu(j)
    175                lonfi((j-2)*iim+1+i)= rlonv(i)
    176             ENDDO
    177          ENDDO
    178          latfi(ngridmx)= rlatu(jjp1)
    179          lonfi(ngridmx)= 0.
    180          
    181          ! build airefi(), mesh area on physics grid
    182          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
    183          ! Poles are single points on physics grid
    184          airefi(1)=airefi(1)*iim
    185          airefi(ngridmx)=airefi(ngridmx)*iim
    186 
    187          CALL inifis(ngridmx,llm,nq,day_ini,daysec,dtphys,
    188      .                latfi,lonfi,airefi,rad,g,r,cpp)
    189       ENDIF
     169!      IF (firstcal) THEN
     170!         latfi(1)=rlatu(1)
     171!         lonfi(1)=0.
     172!         DO j=2,jjm
     173!            DO i=1,iim
     174!               latfi((j-2)*iim+1+i)= rlatu(j)
     175!               lonfi((j-2)*iim+1+i)= rlonv(i)
     176!            ENDDO
     177!         ENDDO
     178!         latfi(ngridmx)= rlatu(jjp1)
     179!         lonfi(ngridmx)= 0.
     180!         
     181!         ! build airefi(), mesh area on physics grid
     182!         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
     183!         ! Poles are single points on physics grid
     184!         airefi(1)=airefi(1)*iim
     185!         airefi(ngridmx)=airefi(ngridmx)*iim
     186!
     187!         CALL inifis(ngridmx,llm,nq,day_ini,daysec,dtphys,
     188!     .                latfi,lonfi,airefi,rad,g,r,cpp)
     189!      ENDIF
    190190
    191191c

  • trunk/LMDZ.MARS/libf/dyn3d/defrun_new.F

    r999 r1130  
    3636c   --------------
    3737! to use  'getin'
    38       USE ioipsl_getincom
     38      use ioipsl_getincom, only: getin
     39      use control_mod, only: ndynstep, day_step, iperiod, iconser,
     40     &                       idissip, iphysiq, anneeref, ecritphy,
     41     &                       ecritstart, timestart, nday_r
    3942      IMPLICIT NONE
    4043
    4144#include "dimensions.h"
    4245#include "paramet.h"
    43 #include "control.h"
     46!#include "control.h"
    4447#include "logic.h"
    4548#include "serre.h"
    … …  
    8386!le modele martien et ne sont donc plus lues dans "run.def"
    8487
    85         anneeref=0
    86         ! Note: anneref is a common in 'control.h'
     88      anneeref=0
    8789
    8890      OPEN(tapedef,file='run.def',status='old',form='formatted'

  • trunk/LMDZ.MARS/libf/dyn3d/dynetat0.F

    r1036 r1130  
    33     
    44      use netcdf
    5       use infotrac, only: tnom
     5      use infotrac, only: tname
     6      use control_mod, only: timestart
    67     
    78      IMPLICIT NONE
    … …  
    3940#include "logic.h"
    4041!#include "advtrac.h"
    41 #include "control.h"
     42!#include "control.h"
    4243
    4344c   Arguments:
    … …  
    381382!           WRITE(str3(2:3),'(i2.2)') iq
    382383!           ierr =  NF_INQ_VARID (nid, str3, nvarid)
    383 ! NB: tracers are now read in using their name ('tnom' from infotrac)
    384 !           write(*,*) "  loading tracer:",trim(tnom(iq))
    385            ierr=nf90_inq_varid(nid,tnom(iq),nvarid)
     384! NB: tracers are now read in using their name ('tname' from infotrac)
     385!           write(*,*) "  loading tracer:",trim(tname(iq))
     386           ierr=nf90_inq_varid(nid,tname(iq),nvarid)
    386387           IF (ierr .NE. nf90_noerr) THEN
    387388!              PRINT*, "dynetat0: Le champ <"//str3//"> est absent"
    388               PRINT*, "dynetat0: Le champ <"//trim(tnom(iq))//
     389              PRINT*, "dynetat0: Le champ <"//trim(tname(iq))//
    389390     &                "> est absent"
    390391              PRINT*, "          Il est donc initialise a zero"
    … …  
    396397             IF (ierr .NE. nf90_noerr) THEN
    397398!                 PRINT*, "dynetat0: Lecture echouee pour "//str3
    398                PRINT*, "dynetat0: Lecture echouee pour "//trim(tnom(iq))
     399               PRINT*,"dynetat0: Lecture echouee pour "//trim(tname(iq))
    399400               CALL abort
    400401             ENDIF

  • trunk/LMDZ.MARS/libf/dyn3d/dynredem.F

    r1106 r1130  
    11      SUBROUTINE dynredem0(fichnom,idayref,anneeref,phis,nq)
    2       use infotrac, only: tnom
     2      use infotrac, only: tname
    33      IMPLICIT NONE
    44c=======================================================================
    … …  
    912912!               ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 12,
    913913!     .                          "Traceurs "//str3)
    914              txt="Traceur "//trim(tnom(iq))
    915 #ifdef NC_DOUBLE
    916                ierr=NF_DEF_VAR(nid,tnom(iq),NF_DOUBLE,4,dims4,nvarid)
    917 #else
    918                ierr=NF_DEF_VAR(nid,tnom(iq),NF_FLOAT,4,dims4,nvarid)
     914             txt="Traceur "//trim(tname(iq))
     915#ifdef NC_DOUBLE
     916               ierr=NF_DEF_VAR(nid,tname(iq),NF_DOUBLE,4,dims4,nvarid)
     917#else
     918               ierr=NF_DEF_VAR(nid,tname(iq),NF_FLOAT,4,dims4,nvarid)
    919919#endif
    920920               ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",
    … …  
    964964      SUBROUTINE dynredem1(fichnom,time,
    965965     .                     vcov,ucov,teta,q,nq,masse,ps)
    966       use infotrac, only: nqtot, tnom
     966      use infotrac, only: nqtot, tname
    967967      IMPLICIT NONE
    968968c=================================================================
    … …  
    11051105!            WRITE(str3(2:3),'(i2.2)') iq
    11061106!            ierr = NF_INQ_VARID(nid, str3, nvarid)
    1107             ierr=NF_INQ_VARID(nid,tnom(iq),nvarid)
     1107            ierr=NF_INQ_VARID(nid,tname(iq),nvarid)
    11081108            IF (ierr .NE. NF_NOERR) THEN
    11091109!               PRINT*, "Variable "//str3//" n est pas definie"
    1110               PRINT*, "Variable "//trim(tnom(iq))//" n est pas definie"
     1110              PRINT*, "Variable "//trim(tname(iq))//" n est pas definie"
    11111111              CALL abort
    11121112            ENDIF

  • trunk/LMDZ.MARS/libf/dyn3d/gcm.F

    r1036 r1130  
    22
    33      use infotrac, only: iniadvtrac, nqtot, iadv
     4      use control_mod, only: day_step, iperiod, iphysiq, ndynstep,
     5     &                       nday_r, idissip, iconser, ecritstart,
     6     &                       ecritphy
     7      use comgeomphy, only: initcomgeomphy
    48      IMPLICIT NONE
    59
    … …  
    4246#include "logic.h"
    4347#include "temps.h"
    44 #include "control.h"
     48!#include "control.h"
    4549#include "ener.h"
    4650#include "netcdf.inc"
    … …  
    135139      logical callgroupeun
    136140      parameter (callgroupeun = .false.)
     141
     142c-----------------------------------------------------------------------
     143c    variables pour l'initialisation de la physique :
     144c    ------------------------------------------------
     145      INTEGER ngridmx
     146      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
     147      REAL zcufi(ngridmx),zcvfi(ngridmx)
     148      REAL latfi(ngridmx),lonfi(ngridmx)
     149      REAL airefi(ngridmx)
     150      SAVE latfi, lonfi, airefi
     151      INTEGER i,j
     152
    137153c-----------------------------------------------------------------------
    138154c   Initialisations:
    … …  
    145161c-----------------------------------------------------------------------
    146162      CALL defrun_new( 99, .TRUE. )
     163
     164!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     165! FH 2008/05/02
     166! A nettoyer. On ne veut qu'une ou deux routines d'interface
     167! dynamique -> physique pour l'initialisation
     168!#ifdef CPP_PHYS
     169      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
     170      call initcomgeomphy
     171!#endif
     172!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    147173
    148174! Initialize tracers
    … …  
    196222c
    197223
     224c-----------------------------------------------------------------------
     225c   Initialisation de la physique :
     226c   -------------------------------
     227
     228!      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
     229         latfi(1)=rlatu(1)
     230         lonfi(1)=0.
     231         zcufi(1) = cu(1)
     232         zcvfi(1) = cv(1)
     233         DO j=2,jjm
     234            DO i=1,iim
     235               latfi((j-2)*iim+1+i)= rlatu(j)
     236               lonfi((j-2)*iim+1+i)= rlonv(i)
     237               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
     238               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
     239            ENDDO
     240         ENDDO
     241         latfi(ngridmx)= rlatu(jjp1)
     242         lonfi(ngridmx)= 0.
     243         zcufi(ngridmx) = cu(ip1jm+1)
     244         zcvfi(ngridmx) = cv(ip1jm-iim)
     245
     246         ! build airefi(), mesh area on physics grid
     247         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
     248         ! Poles are single points on physics grid
     249         airefi(1)=airefi(1)*iim
     250         airefi(ngridmx)=airefi(ngridmx)*iim
     251
     252! Initialisation de la physique: pose probleme quand on tourne
     253! SANS physique, car iniphysiq.F est dans le repertoire phy[]...
     254! Il faut une cle CPP_PHYS
     255!#ifdef CPP_PHYS
     256!         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
     257         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys,
     258     &                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
     259     &                1)
     260!     &                iflag_phys)
     261!#endif
     262!         call_iniphys=.false.
     263!      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
     264!        call inifis(ngridmx,llm,nqtot,day_ini,daysec,dtphys,
     265!     .                latfi,lonfi,airefi,rad,g,r,cpp)
    198266
    199267      call dump2d(iip1,jjp1,ps,'PRESSION SURFACE')

  • trunk/LMDZ.MARS/libf/dyn3d/infotrac.F90

    r1036 r1130  
    55  INTEGER, SAVE :: nqtot
    66  INTEGER,allocatable :: iadv(:)   ! tracer advection scheme number
    7   CHARACTER(len=20),allocatable ::  tnom(:) ! tracer name
     7  CHARACTER(len=20),allocatable ::  tname(:) ! tracer name
    88
    99CONTAINS
    … …  
    1818      IMPLICIT NONE
    1919
    20 #include "dimensions.h"
     20!#include "dimensions.h"
    2121!#include "advtrac.h"
    22 #include "control.h"
     22!#include "control.h"
    2323
    2424! routine arguments:
    … …  
    4747        ! allocate arrays:
    4848        allocate(iadv(nq))
    49         allocate(tnom(nq))
     49        allocate(tname(nq))
    5050       
    5151        ! initialize advection schemes to Van-Leer for all tracers
    … …  
    5656        do iq=1,nq
    5757        ! minimal version, just read in the tracer names, 1 per line
    58           read(90,*,iostat=ierr) tnom(iq)
     58          read(90,*,iostat=ierr) tname(iq)
    5959          if (ierr.ne.0) then
    6060            write(*,*) 'iniadvtrac: error reading tracer names...'

  • trunk/LMDZ.MARS/libf/dyn3d/ini_archive.F

    r1047 r1130  
    4848#include "description.h"
    4949#include "serre.h"
    50 #include "control.h"
     50!#include "control.h"
    5151!#include"comsoil.h"
    5252

  • trunk/LMDZ.MARS/libf/dyn3d/iniconst.F

    r38 r1130  
    11      SUBROUTINE iniconst
    22
     3      use control_mod, only: iphysiq, idissip
    34      IMPLICIT NONE
    45c
    … …  
    1314#include "comconst.h"
    1415#include "temps.h"
    15 #include "control.h"
     16!#include "control.h"
    1617#include "comvert.h"
    1718

  • trunk/LMDZ.MARS/libf/dyn3d/inidissip.F

    r758 r1130  
    88c   -------------
    99
     10      use control_mod, only: idissip, iperiod
    1011      IMPLICIT NONE
    1112#include "dimensions.h"
    … …  
    1415#include "comconst.h"
    1516#include "comvert.h"
    16 #include "control.h"
     17!#include "control.h"
    1718
    1819      LOGICAL lstardis

  • trunk/LMDZ.MARS/libf/dyn3d/lect_start_archive.F

    r1047 r1130  
    1717c
    1818c=======================================================================
    19       use infotrac, only: tnom
     19      use infotrac, only: tname
    2020      use comsoil_h, only: nsoilmx, layer, mlayer, volcapa, inertiedat
    2121      implicit none
    … …  
    3232#include "comvert.h"
    3333#include "comgeom2.h"
    34 #include "control.h"
     34!#include "control.h"
    3535#include "logic.h"
    3636#include "description.h"
    3737#include "ener.h"
    3838#include "temps.h"
    39 #include "lmdstd.h"
     39!#include "lmdstd.h"
    4040#include "netcdf.inc"
    4141!#include "tracer.h"
    … …  
    760760          write(txt,'(a5,i2.2)')'qsurf',iq
    761761        ELSE
    762           txt=trim(tnom(iq))//"_surf"
     762          txt=trim(tname(iq))//"_surf"
    763763          if (txt.eq."h2o_vap") then
    764764            ! There is no surface tracer for h2o_vap;
    … …  
    948948          write(txt,'(a1,i2.2)')'q',iq
    949949        ELSE
    950           txt=tnom(iq)
     950          txt=tname(iq)
    951951        ENDIF
    952952        write(*,*)"lect_start_archive: loading tracer ",trim(txt)

  • trunk/LMDZ.MARS/libf/dyn3d/write_archive.F

    r1047 r1130  
    3838#include "dimphys.h"
    3939#include "paramet.h"
    40 #include "control.h"
     40!#include "control.h"
    4141#include "comvert.h"
    4242#include "comgeom.h"

  • trunk/LMDZ.MARS/libf/dyn3d/writediagdyn.F90

    r410 r1130  
    1212! NB: the rate a which outputs are made can be changed (see parameter isample)
    1313!
     14use control_mod, only: iphysiq, day_step
    1415implicit none
    1516
    1617#include"dimensions.h"
    1718#include"paramet.h"
    18 #include"control.h"
     19!#include"control.h"
    1920#include"netcdf.inc"
    2021

  • trunk/LMDZ.MARS/libf/phymars/albedocaps.F90

    r1047 r1130  
    66! to use the 'getin' routine
    77use ioipsl_getincom, only: getin
    8 #ifdef MESOSCALE
    9 use comgeomfi_h
    10 #endif
     8use comgeomfi_h, only: lati ! grid point latitudes (rad)
    119use surfdat_h, only: TESicealbedo, TESice_Ncoef, TESice_Scoef, &
    1210                     emisice, albedice, watercaptag, albedo_h2o_ice, &
    … …  
    5957
    6058do ig=1,ngrid
    61 #ifndef MESOSCALE
    62   if (ig.gt.ngrid/2+1) then
    63 #else
    64   if (lati(ig)*180./acos(-1.).lt.0.) then
    65 #endif
     59  if (lati(ig).lt.0.) then
    6660    icap=2 ! Southern hemisphere
    6761  else
    … …  
    7468    ! set the surface albedo to be the ice albedo
    7569    if (TESicealbedo) then
    76 !      write(*,*) "albedocaps: call TES_icecap_albedo"
    77 !      write(*,*) "albedocaps: zls=",zls," ig=",ig
    7870      call TES_icecap_albedo(zls,ig,psolaralb(ig,1),icap)
    79 !      write(*,*) "albedocaps: psolaralb(ig,1)=",psolaralb(ig,1)
    8071      psolaralb(ig,2)=psolaralb(ig,1)
    8172    else
    … …  
    8677    ! there is a water ice cap: set the surface albedo to the water ice one
    8778    ! to do : emissivity
    88     !write(*,*) "watercaptag in albedocaps:",ig
    8979    emisref(ig) = 1
    9080    psolaralb(ig,1)=albedo_h2o_ice
    … …  
    10595use comgeomfi_h, only: lati, long
    10696use surfdat_h, only: albedice, TESice_Ncoef, TESice_Scoef
     97use netcdf, only: nf90_open, NF90_NOWRITE, NF90_NOERR, &
     98                  nf90_strerror, nf90_inq_varid, nf90_get_var, nf90_close
     99                 
    107100implicit none
    108101#include"dimensions.h"
    … …  
    110103!#include"surfdat.h"
    111104!#include"comgeomfi.h"
    112 #include"netcdf.inc"
    113105#include"datafile.h"
    114106
    … …  
    169161! Load TES albedoes for Northern Hemisphere
    170162  ! Note: datafile() is defined in "datafile.h"
    171   ierr=NF_OPEN(trim(datafile)//"/npsc_albedo.nc",NF_NOWRITE,nid)
    172   IF (ierr.NE.NF_NOERR) THEN
     163  ierr=nf90_open(trim(datafile)//"/npsc_albedo.nc",NF90_NOWRITE,nid)
     164  IF (ierr.NE.NF90_NOERR) THEN
    173165    write(*,*)'Problem opening npsc_albedo.nc (phymars/albedocaps.F90)'
    174166    write(*,*)'It should be in :',trim(datafile),'/'
    … …  
    179171    write(*,*)'   http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
    180172    CALL ABORT
     173  ELSE
     174    write(*,*) "albedocaps: using file ",trim(datafile)//"/npsc_albedo.nc"
    181175  ENDIF
    182176 
    183   ierr=NF_INQ_VARID(nid,"longitude",nvarid)
    184   if (ierr.ne.NF_NOERR) then
     177  ierr=nf90_inq_varid(nid,"longitude",nvarid)
     178  if (ierr.ne.NF90_NOERR) then
    185179    write(*,*) "Failed to find longitude in file!"
    186   else
    187     ierr=NF_GET_VAR_REAL(nid,nvarid,TESlon)
    188   endif
    189  
    190   ierr=NF_INQ_VARID(nid,"latitude",nvarid)
    191   if (ierr.ne.NF_NOERR) then
     180    write(*,*)trim(nf90_strerror(ierr))
     181    stop
     182  else
     183    ierr=nf90_get_var(nid,nvarid,TESlon)
     184    if (ierr.ne.NF90_NOERR) then
     185      write(*,*) "Failed loading longitude data from file!"
     186      write(*,*)trim(nf90_strerror(ierr))
     187      stop
     188    endif
     189  endif
     190 
     191  ierr=nf90_inq_varid(nid,"latitude",nvarid)
     192  if (ierr.ne.NF90_NOERR) then
    192193    write(*,*) "Failed to find latitude in file!"
    193   else
    194     ierr=NF_GET_VAR_REAL(nid,nvarid,TESlatn)
    195   endif
    196  
    197   ierr=NF_INQ_VARID(nid,"time",nvarid)
    198   if (ierr.ne.NF_NOERR) then
     194    write(*,*)trim(nf90_strerror(ierr))
     195    stop
     196  else
     197    ierr=nf90_get_var(nid,nvarid,TESlatn)
     198    if (ierr.ne.NF90_NOERR) then
     199      write(*,*) "Failed loading latitude data from file!"
     200      write(*,*)trim(nf90_strerror(ierr))
     201      stop
     202    endif
     203  endif
     204 
     205  ierr=nf90_inq_varid(nid,"time",nvarid)
     206  if (ierr.ne.NF90_NOERR) then
    199207    write(*,*) "Failed to find time in file!"
    200   else
    201     ierr=NF_GET_VAR_REAL(nid,nvarid,TESls)
    202   endif
    203  
    204   ierr=NF_INQ_VARID(nid,"albedo",nvarid)
    205   if (ierr.ne.NF_NOERR) then
     208    write(*,*)trim(nf90_strerror(ierr))
     209    stop
     210  else
     211    ierr=nf90_get_var(nid,nvarid,TESls)
     212    if (ierr.ne.NF90_NOERR) then
     213      write(*,*) "Failed loading time data from file!"
     214      write(*,*)trim(nf90_strerror(ierr))
     215      stop
     216    endif
     217  endif
     218 
     219  ierr=nf90_inq_varid(nid,"albedo",nvarid)
     220  if (ierr.ne.NF90_NOERR) then
    206221    write(*,*) "Failed to find albedo in file!"
    207   else
    208     ierr=NF_GET_VAR_REAL(nid,nvarid,TESalbn)
    209   endif
     222    write(*,*)trim(nf90_strerror(ierr))
     223    stop
     224  else
     225    ierr=nf90_get_var(nid,nvarid,TESalbn)
     226    if (ierr.ne.NF90_NOERR) then
     227      write(*,*) "Failed loading albedo data from file!"
     228      write(*,*)trim(nf90_strerror(ierr))
     229      stop
     230    endif
     231  endif
     232
     233  ierr=nf90_close(nid)
    210234
    211235! Load albedoes for Southern Hemisphere
    212   ierr=NF_OPEN(trim(datafile)//"/spsc_albedo.nc",NF_NOWRITE,nid)
    213   IF (ierr.NE.NF_NOERR) THEN
     236  ierr=nf90_open(trim(datafile)//"/spsc_albedo.nc",NF90_NOWRITE,nid)
     237  IF (ierr.NE.NF90_NOERR) THEN
    214238    write(*,*)'Problem opening spsc_albedo.nc (phymars/albedocaps.F90)'
    215239    write(*,*)'It should be in :',trim(datafile),'/'
    … …  
    220244    write(*,*)'   http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
    221245    CALL ABORT
     246  ELSE
     247    write(*,*) "albedocaps: using file ",trim(datafile)//"/spsc_albedo.nc"
    222248  ENDIF
    223249
    224   ierr=NF_INQ_VARID(nid,"latitude",nvarid)
    225   if (ierr.ne.NF_NOERR) then
     250  ierr=nf90_inq_varid(nid,"latitude",nvarid)
     251  if (ierr.ne.NF90_NOERR) then
    226252    write(*,*) "Failed to find latitude in file!"
    227   else
    228     ierr=NF_GET_VAR_REAL(nid,nvarid,TESlats)
    229   endif
    230 
    231   ierr=NF_INQ_VARID(nid,"albedo",nvarid)
    232   if (ierr.ne.NF_NOERR) then
     253    write(*,*)trim(nf90_strerror(ierr))
     254    stop
     255  else
     256    ierr=nf90_get_var(nid,nvarid,TESlats)
     257    if (ierr.ne.NF90_NOERR) then
     258      write(*,*) "Failed loading latitude data from file!"
     259      write(*,*)trim(nf90_strerror(ierr))
     260      stop
     261    endif
     262  endif
     263
     264  ierr=nf90_inq_varid(nid,"albedo",nvarid)
     265  if (ierr.ne.NF90_NOERR) then
    233266    write(*,*) "Failed to find albedo in file!"
    234   else
    235     ierr=NF_GET_VAR_REAL(nid,nvarid,TESalbs)
    236   endif
     267    write(*,*)trim(nf90_strerror(ierr))
     268    stop
     269  else
     270    ierr=nf90_get_var(nid,nvarid,TESalbs)
     271    if (ierr.ne.NF90_NOERR) then
     272      write(*,*) "Failed loading albedo data from file!"
     273      write(*,*)trim(nf90_strerror(ierr))
     274      stop
     275    endif
     276  endif
     277
     278  ierr=nf90_close(nid)
    237279
    238280! constants:

  • trunk/LMDZ.MARS/libf/phymars/co2snow.F

    r1047 r1130  
    33
    44      use surfdat_h, only: iceradius, dtemisice
     5      use comgeomfi_h, only: lati ! grid point latitudes (rad)
    56      IMPLICIT NONE
    67
    … …  
    107108c   dtemis: Time scale for increasing the ice emissivity
    108109
    109            IF(ig.GT.ngrid/2+1) THEN
    110               icap=2
     110           IF(lati(ig).LT. 0.) THEN
     111              icap=2 ! Southern hemisphere
    111112           ELSE
    112               icap=1
     113              icap=1 ! Northern Hemisphere
    113114           ENDIF
    114115

  • trunk/LMDZ.MARS/libf/phymars/comgeomfi_h.F90

    r1047 r1130  
    55
    66       ! These arrays are allocated in inifis
    7        REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: long,lati,area
     7       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: long ! longitudes (rad)
     8       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: lati ! latitudes (rad)
     9       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: area ! mesh area (m2)
    810       REAL,SAVE :: totarea
    911

  • trunk/LMDZ.MARS/libf/phymars/datareadnc.F

    r690 r1130  
    9999
    100100
    101 #include "lmdstd.h"
     101!#include "lmdstd.h"
    102102#include "fxyprim.h"
    103103

  • trunk/LMDZ.MARS/libf/phymars/eofdump_mod.F90

    r1047 r1130  
    2727      integer,intent(in) :: ngrid ! total number of physics grid points
    2828      integer,intent(in) :: nlayer ! number of atmospheric layers
    29       real*4 u(ngrid,nlayer)
    30       real*4 v(ngrid,nlayer)
    31       real*4 t(ngrid,nlayer)
    32       real*4 rho(ngrid,nlayer)
    33       real*4 ps(ngrid)
     29      real,intent(in) :: u(ngrid,nlayer)
     30      real,intent(in) :: v(ngrid,nlayer)
     31      real,intent(in) :: t(ngrid,nlayer)
     32      real,intent(in) :: rho(ngrid,nlayer)
     33      real,intent(in) :: ps(ngrid)
    3434      integer,save :: count=0
    3535      integer i,j,l, ig
    … …  
    5656          do j=1+eofskip/2,jjm+1,eofskip
    5757            ig = 1+ (j-2)*iim +i
     58#ifdef NC_DOUBLE
     59            write(uedata) (real(u(ig,l)),l=1,nlayer)
     60            write(uedata) (real(v(ig,l)),l=1,nlayer)
     61            write(uedata) (real(t(ig,l)),l=1,nlayer)
     62            write(uedata) (real(rho(ig,l)),l=1,nlayer)
     63            write(uedata) real(ps(ig))
     64#else
    5865            write(uedata) (u(ig,l),l=1,nlayer)
    5966            write(uedata) (v(ig,l),l=1,nlayer)
    … …  
    6168            write(uedata) (rho(ig,l),l=1,nlayer)
    6269            write(uedata) ps(ig)
     70#endif
    6371          enddo
    6472        enddo
    … …  
    112120          if(j.eq.1) stop 'Problem in ineofdump.F'
    113121          if(j.eq.jjm+1) stop 'Problem in ineofdump.F'
     122#ifdef NC_DOUBLE
     123          write(uehead,*) real(long(ig)*180./pi),real(lati(ig)*180./pi)
     124#else
    114125          write(uehead,*) long(ig)*180./pi, lati(ig)*180./pi
     126#endif
    115127!         write(*,*) 'eof grid j=',j,' lat= ', lati(ig)*180./pi
    116128        enddo
    117129      enddo
    118130
     131#ifdef NC_DOUBLE
     132      write(uehead,*) real(aps)
     133      write(uehead,*) real(bps)
     134#else
    119135      write(uehead,*) aps
    120136      write(uehead,*) bps
     137#endif
    121138      close(uehead)
    122139!

  • trunk/LMDZ.MARS/libf/phymars/inifis.F

    r1112 r1130  
    5959      use yomlw_h, only: ini_yomlw_h
    6060      use conc_mod, only: ini_conc_mod
     61      use control_mod, only: ecritphy
    6162
    6263#ifdef MESOSCALE
    … …  
    136137!      ENDIF
    137138
     139      ! read in 'ecritphy' (frequency of calls to physics, in dynamical steps)
     140      ! (also done in dyn3d/defrun_new but not in LMDZ.COMMON)
     141      call getin("ecritphy",ecritphy)
     142     
    138143! --------------------------------------------------------------
    139144!  Reading the "callphys.def" file controlling some key options

  • trunk/LMDZ.MARS/libf/phymars/inistats.F

    r690 r1130  
    11      subroutine inistats(ierr)
    22
     3      use mod_phys_lmdz_para, only : is_master
    34      implicit none
    45
    … …  
    4243         pseudoalt(l)=-10.*log(presnivs(l)/preff)   
    4344      enddo
     45
     46      if (is_master) then
     47      ! only the master needs do this
    4448
    4549      ierr = NF_CREATE("stats.nc",IOR(NF_CLOBBER,NF_64BIT_OFFSET),nid)
    … …  
    110114      ierr=NF_CLOSE(nid)
    111115
     116      endif ! of if (is_master)
    112117      end subroutine inistats
    113118

  • trunk/LMDZ.MARS/libf/phymars/initracer.F

    r1047 r1130  
    22
    33#ifndef MESOSCALE
    4        use infotrac, only: tnom
     4       use infotrac, only: tname
    55#endif
    66       use tracer_mod
    … …  
    8181        txt=" "
    8282        write(txt,'(a1,i2.2)') 'q',iq
    83         if (txt.eq.tnom(iq)) then
     83        if (txt.eq.tname(iq)) then
    8484          count=count+1
    8585        endif
    … …  
    9595      ! copy tracer names from dynamics
    9696      do iq=1,nq
    97         noms(iq)=tnom(iq)
     97        noms(iq)=tname(iq)
    9898      enddo
    9999#endif

  • trunk/LMDZ.MARS/libf/phymars/mkstat.F90

    r38 r1130  
    1010!  Yann W. july 2003
    1111
     12use mod_phys_lmdz_para, only : is_master
    1213
    1314implicit none
    … …  
    3637! Incrementation of count for the last step, which is not done in wstats
    3738count(istime)=count(istime)+1
     39
     40if (is_master) then
     41! only the master needs do this
    3842
    3943ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
    … …  
    162166ierr= NF_CLOSE(nid)
    163167
     168endif ! of if (is_master)
     169
    164170end

  • trunk/LMDZ.MARS/libf/phymars/newcondens.F

    r1114 r1130  
    217217      ENDIF ! of IF (firstcall)
    218218      zcpi=1./cpp
     219
    219220c
    220221c======================================================================

  • trunk/LMDZ.MARS/libf/phymars/newstart.F

    r1087 r1130  
    1717! to use  'getin'
    1818      use ioipsl_getincom, only: getin
    19       use infotrac, only: iniadvtrac, nqtot, tnom
    20       use tracer_mod, only: noms, igcm_h2o_vap, igcm_h2o_ice
     19      use infotrac, only: iniadvtrac, nqtot, tname
     20      use tracer_mod, only: noms, igcm_dust_number, igcm_dust_mass,
     21     &                      igcm_ccn_number, igcm_ccn_mass,
     22     &                      igcm_h2o_vap, igcm_h2o_ice
    2123      use surfdat_h, only: phisfi, z0, zmea, zstd, zsig, zgam, zthe,
    2224     &                     albedodat, z0_default
    2325      use comsoil_h, only: inertiedat, layer, mlayer, nsoilmx
     26      use control_mod, only: day_step, iphysiq, anneeref
     27      use phyredem, only: physdem0, physdem1
     28      use iostart, only: open_startphy
     29      use comgeomphy, only: initcomgeomphy
    2430      implicit none
    2531
    … …  
    3642#include "comvert.h"
    3743#include "comgeom2.h"
    38 #include "control.h"
     44!#include "control.h"
    3945#include "logic.h"
    4046#include "description.h"
    4147#include "ener.h"
    4248#include "temps.h"
    43 #include "lmdstd.h"
     49!#include "lmdstd.h"
    4450#include "comdissnew.h"
    4551#include "clesph0.h"
    … …  
    142148c variables diverses
    143149c-------------------
    144       real choix_1
     150      real choix_1 ! ==0 : read start_archive file ; ==1: read start files
    145151      character*80      fichnom
    146152      integer Lmodif,iq
    … …  
    183189      preff  = 610.    ! for Mars, instead of 101325. (Earth)
    184190      pa= 20           ! for Mars, instead of 500 (Earth)
     191
     192! initialize "serial/parallel" related stuff
     193      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
     194      call initcomgeomphy
    185195
    186196! Load tracer number and names:
    … …  
    316326c-----------------------------------------------------------------------
    317327      if (choix_1.eq.0) then
     328         ! tabfi requires that input file be first opened by open_startphy(fichnom)
     329         fichnom = 'start_archive.nc'
     330         call open_startphy(fichnom)
    318331         call tabfi (nid,Lmodif,tab0,day_ini,lllm,p_rad,
    319332     .            p_omeg,p_g,p_mugaz,p_daysec,time)
    320333      else if (choix_1.eq.1) then
     334         fichnom = 'startfi.nc'
     335         call open_startphy(fichnom)
    321336         call tabfi (nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
    322337     .            p_omeg,p_g,p_mugaz,p_daysec,time)
    … …  
    461476        txt=" "
    462477        write(txt,'(a1,i2.2)') 'q',iq
    463         if (txt.eq.tnom(iq)) then
     478        if (txt.eq.tname(iq)) then
    464479          count=count+1
    465480        endif
    … …  
    471486      if (count.eq.nqtot) then
    472487        write(*,*) 'Newstart: updating tracer names'
    473         ! copy noms(:) to tnom(:) to have matching tracer names in physics
     488        ! copy noms(:) to tname(:) to have matching tracer names in physics
    474489        ! and dynamics
    475         tnom(1:nqtot)=noms(1:nqtot)
     490        tname(1:nqtot)=noms(1:nqtot)
    476491      endif
    477492
    … …  
    497512     $ tracer'
    498513      write(*,*) 'q=profile    : specify a profile for a tracer'
     514      write(*,*) 'freedust     : rescale dust to a true value'
    499515      write(*,*) 'ini_q        : tracers initialization for chemistry
    500516     $ and water vapour'
    … …  
    723739          write(*,*) 'Which tracer name do you want to change ?'
    724740          do iq=1,nqtot
    725             write(*,'(i3,a3,a20)')iq,' : ',trim(tnom(iq))
     741            write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
    726742          enddo
    727743          write(*,'(a35,i3)')
    … …  
    730746          read(*,*) iq
    731747          if ((iq.ge.1).and.(iq.le.nqtot)) then
    732             write(*,*)'Change tracer name ',trim(tnom(iq)),' to ?'
     748            write(*,*)'Change tracer name ',trim(tname(iq)),' to ?'
    733749            read(*,*) txt
    734             tnom(iq)=txt
     750            tname(iq)=txt
    735751            write(*,*)'Do you want to change another tracer name (y/n)?'
    736752            read(*,'(a)') yes
    … …  
    768784             write(*,*) 'Which tracer do you want to modify ?'
    769785             do iq=1,nqtot
    770                write(*,*)iq,' : ',trim(tnom(iq))
     786               write(*,*)iq,' : ',trim(tname(iq))
    771787             enddo
    772788             write(*,*) '(choose between 1 and ',nqtot,')'
    … …  
    777793               cycle
    778794             endif
    779              write(*,*)'mixing ratio of tracer ',trim(tnom(iq)),
     795             write(*,*)'mixing ratio of tracer ',trim(tname(iq)),
    780796     &                 ' ? (kg/kg)'
    781797             read(*,*) val
    … …  
    787803               ENDDO
    788804             ENDDO
    789              write(*,*) 'SURFACE value of tracer ',trim(tnom(iq)),
     805             write(*,*) 'SURFACE value of tracer ',trim(tname(iq)),
    790806     &                   ' ? (kg/m2)'
    791807             read(*,*) val
    … …  
    804820             write(*,*) 'Which tracer do you want to set?'
    805821             do iq=1,nqtot
    806                write(*,*)iq,' : ',trim(tnom(iq))
     822               write(*,*)iq,' : ',trim(tname(iq))
    807823             enddo
    808824             write(*,*) '(choose between 1 and ',nqtot,')'
    … …  
    814830             endif
    815831             ! look for input file 'profile_tracer'
    816              txt="profile_"//trim(tnom(iq))
     832             txt="profile_"//trim(tname(iq))
    817833             open(41,file=trim(txt),status='old',form='formatted',
    818834     &            iostat=ierr)
    … …  
    831847                   q(:,:,l,iq)=profile(l+1)
    832848                 enddo
    833                  write(*,*)'OK, tracer ',trim(tnom(iq)),' initialized ',
    834      &                     'using values from file ',trim(txt)
     849                 write(*,*)'OK, tracer ',trim(tname(iq)),
     850     &               ' initialized ','using values from file ',trim(txt)
    835851               else
    836852                 write(*,*)'problem reading file ',trim(txt),' !'
    837                  write(*,*)'No modifications to tracer ',trim(tnom(iq))
     853                 write(*,*)'No modifications to tracer ',trim(tname(iq))
    838854               endif
    839855             else
    840856               write(*,*)'Could not find file ',trim(txt),' !'
    841                write(*,*)'No modifications to tracer ',trim(tnom(iq))
     857               write(*,*)'No modifications to tracer ',trim(tname(iq))
    842858             endif
    843859             
     860c       q=profile : initialize tracer with a given profile
     861c       --------------------------------------------------
     862        else if (trim(modif) .eq. 'freedust') then
     863          do l=1,llm
     864            do j=1,jjp1
     865              do i=1,iip1
     866                if (igcm_dust_number .ne. 0)
     867     &          q(i,j,l,igcm_dust_number)=
     868     &                q(i,j,l,igcm_dust_number) * 1e-3 ! grosso modo
     869                if (igcm_dust_mass .ne. 0)
     870     &          q(i,j,l,igcm_dust_mass)=
     871     &                q(i,j,l,igcm_dust_mass)   * 1e-3 ! grosso modo
     872                if (igcm_ccn_number .ne. 0)
     873     &          q(i,j,l,igcm_ccn_number)=
     874     &                q(i,j,l,igcm_ccn_number) * 1e-3 ! grosso modo
     875                if (igcm_ccn_mass .ne. 0)
     876     &          q(i,j,l,igcm_ccn_mass)=
     877     &                q(i,j,l,igcm_ccn_mass)   * 1e-3 ! grosso modo
     878              end do
     879            end do
     880          end do
     881
     882         ! We want to have the very same value at lon -180 and lon 180
     883          do l = 1,llm
     884             do j = 1,jjp1
     885                do iq = 1,nqtot
     886                   q(iip1,j,l,iq) = q(1,j,l,iq)
     887                end do
     888             end do
     889          end do
    844890
    845891c       ini_q : Initialize tracers for chemistry

  • trunk/LMDZ.MARS/libf/phymars/physiq.F

    r1124 r1130  
    2828      use slope_mod, only: theta_sl, psi_sl
    2929      use conc_mod, only: rnew, cpnew, mmean
     30      use control_mod, only: iphysiq, day_step, ecritstart
     31      use phyredem, only: physdem0, physdem1
    3032
    3133#ifdef MESOSCALE
    … …  
    144146#include "planete.h"
    145147!#include "comsaison.h"
    146 #include "control.h"
     148!#include "control.h"
    147149!#include "dimradmars.h"
    148150! naerkind is set in scatterers.h (built when compiling with makegcm -s #)
    … …  
    515517#ifndef MESOSCALE
    516518         if (callslope) call getslopes(ngrid,phisfi)
    517                            
     519
     520         if (ngrid.ne.1) then ! no need to create a restart file in 1d
    518521         call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq,
    519522     .                 ptimestep,pday,time_phys,area,
    520523     .                 albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
    521      
     524         endif
    522525#endif
    523526                   
    … …  
    676679c          Radiative transfer
    677680c          ------------------
    678 
    679681           CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
    680682     $     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
    … …  
    10361038     &                nq,tau,tauscaling,rdust,rice,nuice,
    10371039     &                rsedcloud,rhocloud)
    1038 
     1040     
    10391041c Temperature variation due to latent heat release
    10401042           if (activice) then
    … …  
    11351137c        -------------
    11361138         IF (sedimentation) THEN
    1137            !call zerophys(ngrid*nlayer*nq, zdqsed)
    11381139           zdqsed(1:ngrid,1:nlayer,1:nq)=0
    1139            !call zerophys(ngrid*nq, zdqssed)
    11401140           zdqssed(1:ngrid,1:nq)=0
    11411141
    … …  
    20732073     &                       'Mean reff',
    20742074     &                       'm',2,rave)
    2075             CALL WRITEDIAGFI(ngrid,"Nccntot",
     2075            if (scavenging) then
     2076              CALL WRITEDIAGFI(ngrid,"Nccntot",
    20762077     &                    "condensation nuclei","Nbr/m2",
    20772078     &                    2,Nccntot)
    2078             CALL WRITEDIAGFI(ngrid,"Mccntot",
     2079              CALL WRITEDIAGFI(ngrid,"Mccntot",
    20792080     &                    "mass condensation nuclei","kg/m2",
    20802081     &                    2,Mccntot)
     2082            endif
    20812083            call WRITEDIAGFI(ngrid,'rice','Ice particle size',
    20822084     &                       'm',3,rice)

  • trunk/LMDZ.MARS/libf/phymars/planete.h

    r224 r1130  
    1 c-----------------------------------------------------------------------
    2 c INCLUDE planet.h
     1!-----------------------------------------------------------------------
     2! INCLUDE planet.h
    33
    44      COMMON/planet/aphelie,periheli,year_day,peri_day,                 &
    … …  
    1212     &       timeperi,e_elips,p_elips,unitastr
    1313
    14 c-----------------------------------------------------------------------
     14!-----------------------------------------------------------------------

  • trunk/LMDZ.MARS/libf/phymars/soil_settings.F

    r1047 r1130  
    11      subroutine soil_settings(nid,ngrid,nsoil,tsurf,tsoil,indextime)
    22
    3       use netcdf
     3!      use netcdf
    44      use comsoil_h, only: layer, mlayer, inertiedat, volcapa
     5      use iostart, only: inquire_field_ndims, get_var, get_field,
     6     &                   inquire_field, inquire_dimension_length
    57      implicit none
    68
    … …  
    7476      real,parameter :: default_volcapa=1.e6
    7577
     78      logical :: found,ok
     79     
    7680!======================================================================
    7781
    … …  
    8185! 1.1 Start by reading how many layers of soil there are
    8286
    83         ierr=nf90_inq_dimid(nid,"subsurface_layers",dimid)
    84         if (ierr.ne.nf90_noerr) then
    85          write(*,*)'soil_settings: Problem reading <subsurface_layers>'
    86          write(*,*)trim(nf90_strerror(ierr))
    87          call abort
    88         endif
    89 
    90         ierr=nf90_inquire_dimension(nid,dimid,len=dimlen)
    91         if (ierr.ne.nf90_noerr) then
    92          write(*,*)'soil_settings: Problem getting <subsurface_layers>',
    93      &             'length'
    94          write(*,*)trim(nf90_strerror(ierr))
    95          call abort
    96         endif
     87!       ierr=nf90_inq_dimid(nid,"subsurface_layers",dimid)
     88!       if (ierr.ne.nf90_noerr) then
     89!        write(*,*)'soil_settings: Problem reading <subsurface_layers>'
     90!         write(*,*)trim(nf90_strerror(ierr))
     91!        call abort
     92!       endif
     93
     94!       ierr=nf90_inquire_dimension(nid,dimid,len=dimlen)
     95!       if (ierr.ne.nf90_noerr) then
     96!        write(*,*)'soil_settings: Problem getting <subsurface_layers>',
     97!     &             'length'
     98!        write(*,*)trim(nf90_strerror(ierr))
     99!         call abort
     100!       endif
     101        dimlen=inquire_dimension_length("subsurface_layers")
    97102
    98103        if (dimlen.ne.nsoil) then
    … …  
    107112!     (in ye old days, thermal inertia was only given at the "surface")
    108113      ! Look for thermal inertia data
    109       ierr=nf90_inq_varid(nid,"inertiedat",nvarid)
    110       if (ierr.NE.nf90_noerr) then
    111          write(*,*)'soil_settings: Field <inertiedat> not found!'
    112          write(*,*)trim(nf90_strerror(ierr))
    113          call abort
    114       endif
    115 
    116       ! Read the # of dimensions <inertidat> was defined as using
    117       ierr=nf90_inquire_variable(nid,nvarid,ndims=ndims)
    118       ! if (ndims.eq.1) then we have the "old 2D-surface" format
     114!      ierr=nf90_inq_varid(nid,"inertiedat",nvarid)
     115!      if (ierr.NE.nf90_noerr) then
     116!         write(*,*)'soil_settings: Field <inertiedat> not found!'
     117!         write(*,*)trim(nf90_strerror(ierr))
     118!         call abort
     119!      endif
     120!
     121!      ! Read the # of dimensions <inertidat> was defined as using
     122!      ierr=nf90_inquire_variable(nid,nvarid,ndims=ndims)
     123!      ! if (ndims.eq.1) then we have the "old 2D-surface" format
     124      ndims=inquire_field_ndims("inertiedat")
    119125
    120126! 1.3 Read depths values or set olddepthdef flag and values
    … …  
    136142        enddo
    137143      else ! Look for depth
    138         ierr=nf90_inq_varid(nid,"soildepth",nvarid)
    139         if (ierr.ne.nf90_noerr) then
    140           write(*,*)'soil_settings: Field <soildepth> not found!'
    141           write(*,*)trim(nf90_strerror(ierr))
    142           call abort
    143         endif
     144!        ierr=nf90_inq_varid(nid,"soildepth",nvarid)
     145!        if (ierr.ne.nf90_noerr) then
     146!          write(*,*)'soil_settings: Field <soildepth> not found!'
     147!          write(*,*)trim(nf90_strerror(ierr))
     148!         call abort
     149!        endif
    144150        ! read <depth> coordinate
    145151        if (interpol) then !put values in oldmlayer
    146           ierr=nf90_get_var(nid,nvarid,oldmlayer)
    147           if (ierr.ne.nf90_noerr) then
    148            write(*,*)'soil_settings: Problem while reading <soildepth>'
    149            write(*,*)trim(nf90_strerror(ierr))
    150            call abort
     152!          ierr=nf90_get_var(nid,nvarid,oldmlayer)
     153!          if (ierr.ne.nf90_noerr) then
     154!           write(*,*)'soil_settings: Problem while reading <soildepth>'
     155!           write(*,*)trim(nf90_strerror(ierr))
     156!           call abort
     157!          endif
     158          call get_var("soildepth",oldmlayer,found)
     159          if (.not.found) then
     160            write(*,*)'soil_settings: Problem while reading <soildepth>'
    151161          endif
    152162        else ! put values in mlayer
    153           ierr=nf90_get_var(nid,nvarid,mlayer)
    154           if (ierr.ne.nf90_noerr) then
    155            write(*,*)'soil_settings: Problem while reading <soildepth>'
    156            write(*,*)trim(nf90_strerror(ierr))
    157            call abort
     163!          ierr=nf90_get_var(nid,nvarid,mlayer)
     164!          if (ierr.ne.nf90_noerr) then
     165!           write(*,*)'soil_settings: Problem while reading <soildepth>'
     166!           write(*,*)trim(nf90_strerror(ierr))
     167!           call abort
     168!          endif
     169          call get_var("soildepth",mlayer,found)
     170          if (.not.found) then
     171            write(*,*)'soil_settings: Problem while reading <soildepth>'
    158172          endif
    159173        endif !of if (interpol)
    … …  
    217231
    218232! 3.1 Look (again) for thermal inertia data (to reset nvarid)
    219       ierr=nf90_inq_varid(nid,"inertiedat",nvarid)
    220       if (ierr.NE.nf90_noerr) then
    221          write(*,*)'soil_settings: Field <inertiedat> not found!'
    222          write(*,*)trim(nf90_strerror(ierr))
    223          call abort
    224       endif
     233!      ierr=nf90_inq_varid(nid,"inertiedat",nvarid)
     234!      if (ierr.NE.nf90_noerr) then
     235!         write(*,*)'soil_settings: Field <inertiedat> not found!'
     236!         write(*,*)trim(nf90_strerror(ierr))
     237!         call abort
     238!      endif
    225239
    226240! 3.2 Knowing the # of dimensions <inertidat> was defined as using,
    … …  
    232246       ! Read Surface thermal inertia
    233247       allocate(surfinertia(ngrid))
    234        ierr=nf90_get_var(nid,nvarid,surfinertia)
    235         if (ierr.NE.nf90_noerr) then
    236          write(*,*)'soil_settings: Problem while reading <inertiedat>'
    237          write(*,*)trim(nf90_strerror(ierr))
     248!       ierr=nf90_get_var(nid,nvarid,surfinertia)
     249!        if (ierr.NE.nf90_noerr) then
     250!         write(*,*)'soil_settings: Problem while reading <inertiedat>'
     251!         write(*,*)trim(nf90_strerror(ierr))
     252!         call abort
     253!        endif
     254       call get_field("inertiedat",surfinertia,found)
     255       if (.not.found) then
     256         write(*,*) "soil_settings: Failed loading <inertiedat>"
    238257         call abort
    239         endif
     258       endif
    240259       
    241260       write(*,*)' => Building soil thermal inertia (using reference sur
    … …  
    255274            stop
    256275           endif
    257          endif
    258          ierr=nf90_get_var(nid,nvarid,oldinertiedat)
    259         if (ierr.NE.nf90_noerr) then
    260          write(*,*)'soil_settings: Problem while reading <inertiedat>'
    261          write(*,*)trim(nf90_strerror(ierr))
    262          call abort
     276         endif ! of if (.not.allocated(oldinertiedat))
     277!         ierr=nf90_get_var(nid,nvarid,oldinertiedat)
     278!        if (ierr.NE.nf90_noerr) then
     279!         write(*,*)'soil_settings: Problem while reading <inertiedat>'
     280!         write(*,*)trim(nf90_strerror(ierr))
     281!         call abort
     282!        endif
     283        call get_field("inertiedat",oldinertiedat,found)
     284        if (.not.found) then
     285          write(*,*) "soil_settings: Failed loading <inertiedat>"
     286          call abort
    263287        endif
    264288       else ! put values in therm_i
    265         ierr=nf90_get_var(nid,nvarid,inertiedat)
    266         if (ierr.NE.nf90_noerr) then
    267          write(*,*)'soil_settings: Problem while reading <inertiedat>'
    268          write(*,*)trim(nf90_strerror(ierr))
    269          call abort
    270         endif
     289!        ierr=nf90_get_var(nid,nvarid,inertiedat)
     290!        if (ierr.NE.nf90_noerr) then
     291!         write(*,*)'soil_settings: Problem while reading <inertiedat>'
     292!         write(*,*)trim(nf90_strerror(ierr))
     293!         call abort
     294         call get_field("inertiedat",inertiedat,found)
     295         if (.not.found) then
     296           write(*,*) "soil_settings: Failed loading <inertiedat>"
     297           call abort
     298         endif
     299!        endif
    271300       endif ! of if (interpol)
    272301      endif ! of if (ndims.eq.1)
    … …  
    275304! -------------------------
    276305
    277       ierr=nf90_inq_varid(nid,"tsoil",nvarid)
    278       if (ierr.ne.nf90_noerr) then
     306!      ierr=nf90_inq_varid(nid,"tsoil",nvarid)
     307      ok=inquire_field("tsoil")
     308!      if (ierr.ne.nf90_noerr) then
     309      if (.not.ok) then
    279310        write(*,*)'soil_settings: Field <tsoil> not found!'
    280311        write(*,*)' => Building <tsoil> from surface values <tsurf>'
    … …  
    292323           endif
    293324         endif
    294         ierr=nf90_get_var(nid,nvarid,oldtsoil)
    295         if (ierr.ne.nf90_noerr) then
    296          write(*,*)'soil_settings: Problem while reading <tsoil>'
    297          write(*,*)trim(nf90_strerror(ierr))
    298          call abort
    299         endif
     325!        ierr=nf90_get_var(nid,nvarid,oldtsoil)
     326!        if (ierr.ne.nf90_noerr) then
     327!        write(*,*)'soil_settings: Problem while reading <tsoil>'
     328!         write(*,*)trim(nf90_strerror(ierr))
     329!        call abort
     330!       endif
     331         call get_field("tsoil",oldtsoil,found)
     332         if (.not.found) then
     333           write(*,*) "soil_settings: Failed loading <tsoil>"
     334           call abort
     335         endif
    300336       else ! put values in tsoil
    301         corner(1)=1
    302         corner(2)=1
    303         corner(3)=indextime
    304         edges(1)=ngrid
    305         edges(2)=nsoil
    306         edges(3)=1
    307         !ierr=nf90_get_var(nid,nvarid,tsoil,corner,edges)
    308         ierr=nf90_get_var(nid,nvarid,tsoil)
    309         if (ierr.ne.nf90_noerr) then
    310          write(*,*)'soil_settings: Problem while reading <tsoil>'
    311          write(*,*)trim(nf90_strerror(ierr))
    312          call abort
    313         endif
     337!        corner(1)=1
     338!        corner(2)=1
     339!        corner(3)=indextime
     340!        edges(1)=ngrid
     341!        edges(2)=nsoil
     342!        edges(3)=1
     343!        !ierr=nf90_get_var(nid,nvarid,tsoil,corner,edges)
     344!        ierr=nf90_get_var(nid,nvarid,tsoil)
     345!        if (ierr.ne.nf90_noerr) then
     346!        write(*,*)'soil_settings: Problem while reading <tsoil>'
     347!         write(*,*)trim(nf90_strerror(ierr))
     348!        call abort
     349!       endif
     350         call get_field("tsoil",tsoil,found,timeindex=indextime)
     351         if (.not.found) then
     352           write(*,*) "soil_settings: Failed loading <tsoil>"
     353           call abort
     354         endif
    314355       endif ! of if (interpol)
    315       endif
     356      endif! of if (.not.ok)
    316357
    317358! 5. If necessary, interpolate soil temperatures and thermal inertias

  • trunk/LMDZ.MARS/libf/phymars/start2archive.F

    r1087 r1130  
    1919c=======================================================================
    2020
    21       use infotrac, only: iniadvtrac, nqtot, tnom
     21      use infotrac, only: iniadvtrac, nqtot, tname
    2222      use comsoil_h, only: nsoilmx, inertiedat
    2323      use surfdat_h, only: ini_surfdat_h
    2424      use comsoil_h, only: ini_comsoil_h
     25      use comgeomphy, only: initcomgeomphy
    2526      implicit none
    2627
    … …  
    3435#include "logic.h"
    3536#include "temps.h"
    36 #include "control.h"
     37!#include "control.h"
    3738#include "ener.h"
    3839#include "description.h"
    … …  
    119120      CALL defrun_new(99, .TRUE. )
    120121      grireg   = .TRUE.
     122! initialize "serial/parallel" related stuff
     123      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
     124      call initcomgeomphy
    121125
    122126c=======================================================================
    … …  
    355359c        call write_archive(nid,ntime,'q'//str2,'tracer','kg/kg',
    356360c     .         3,q(1,1,iq))
    357         call write_archive(nid,ntime,tnom(iq),'tracer','kg/kg',
     361        call write_archive(nid,ntime,tname(iq),'tracer','kg/kg',
    358362     &         3,q(1,1,iq))
    359363      end do
    … …  
    365369c       call write_archive(nid,ntime,'qsurf'//str2,'Tracer on surface',
    366370c     $  'kg.m-2',2,qsurfS(1,iq))
    367         txt=trim(tnom(iq))//"_surf"
     371        txt=trim(tname(iq))//"_surf"
    368372        call write_archive(nid,ntime,txt,'Tracer on surface',
    369373     &  'kg.m-2',2,qsurfS(1,iq))
    … …  
    404408c-----------------------------------------------------------------------
    405409
     410      write(*,*) "startarchive: all is well that ends well"
     411     
    406412      end

  • trunk/LMDZ.MARS/libf/phymars/surfini.F

    r1087 r1130  
    88     &                     albedo_h2o_ice, inert_h2o_ice, albedodat,
    99     &                     albedice
     10      use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
     11      use mod_phys_lmdz_para, only : is_master, gather, scatter
    1012      IMPLICIT NONE
    1113c=======================================================================
    … …  
    2729#include "datafile.h"
    2830
    29       INTEGER ngrid,ig,icap,iq,alternate
    30       REAL  piceco2(ngrid),psolaralb(ngrid,2)
    31       REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2)
     31      integer,intent(in) :: ngrid ! number of atmospheric columns
     32      real,intent(in) :: piceco2(ngrid) ! CO2 ice thickness
     33      real,intent(inout) :: qsurf(ngrid,nqmx) ! tracer on surface (kg/m2)
     34      real,intent(out) :: psolaralb(ngrid,2) ! albedo
     35
     36      INTEGER ig,icap,iq,alternate
    3237      REAL icedryness ! ice dryness
    3338     
    3439      ! longwatercaptag is watercaptag. Trick for some compilers
    3540      LOGICAL, DIMENSION(100000) :: longwatercaptag
    36 
    37       EXTERNAL ISMIN,ISMAX
    38       INTEGER ISMIN,ISMAX
    3941     
    4042! There are 3 different modes for ice distribution:
    … …  
    5355      REAL        zelat,zelon
    5456
     57#ifdef CPP_PARA
    5558      INTEGER nb_ice(ngrid,2)              ! number of counts | detected ice for GCM grid
     59#else
     60      INTEGER nb_ice(klon_glo,2)
     61#endif
    5662      INTEGER latice(jjm,2),lonice (iim,2) ! number of counts | detected ice along lat & lon axis
    5763
    … …  
    6470     
    6571      character (len=100) :: zedatafile
     72
     73      ! to handle parallel cases
     74#if CPP_PARA
     75      logical watercaptag_glo(klon_glo)
     76      real dryness_glo(klon_glo)
     77      real lati_glo(klon_glo)
     78      real long_glo(klon_glo)
     79#else
     80      logical watercaptag_glo(ngrid)
     81      real dryness_glo(ngrid)
     82      real lati_glo(ngrid)
     83      real long_glo(ngrid)
     84#endif
     85
    6686c
    6787c=======================================================================
     88! Initialize watercaptag (default is false)
     89      watercaptag_glo(:)=.false.
    6890
    6991c     water ice outliers
    … …  
    93115         endif
    94116         
    95          write(*,*) "Ice dryness ?"
     117         write(*,*) "surfini: Ice dryness ?"
    96118         icedryness=1. ! default value
    97119         call getin("icedryness",icedryness)
    98          write(*,*) " icedryness = ",icedryness
     120         write(*,*) "surfini: icedryness = ",icedryness
    99121         dryness (:) = icedryness
    100122         
    … …  
    125147#else
    126148
    127 
    128 
    129       IF (ngrid .eq. 1) THEN ! special case for 1d --> do nothing
     149      ! To be able to run in parallel, we work on the full grid
     150      ! and dispatch results afterwards
     151
     152      ! start by geting latitudes and logitudes on full grid
     153      ! (in serial mode, this is just a copy)
     154      call gather(lati,lati_glo)
     155      call gather(long,long_glo)
     156
     157      if (is_master) then
     158
     159        IF (ngrid .eq. 1) THEN ! special case for 1d --> do nothing
    130160     
    131161         print*, 'ngrid = 1, do no put ice caps in surfini.F'
    132162
    133       ELSE IF (icelocationmode .eq. 1) THEN
     163        ELSE IF (icelocationmode .eq. 1) THEN
    134164     
    135165         print*,'Surfini: ice caps defined from surface.nc'
    … …  
    144174
    145175           
    146         zedatafile = trim(datafile)
     176         zedatafile = trim(datafile)
    147177 
    148178       
    149         ierr=nf90_open(trim(zedatafile)//'/surface.nc',
    150      &  NF90_NOWRITE,nid)
     179         ierr=nf90_open(trim(zedatafile)//'/surface.nc',
     180     &   NF90_NOWRITE,nid)
    151181     
    152       IF (ierr.NE.nf90_noerr) THEN
     182         IF (ierr.NE.nf90_noerr) THEN
    153183       write(*,*)'Error : cannot open file surface.nc '
    154184       write(*,*)'(in phymars/surfini.F)'
    … …  
    160190       write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
    161191       CALL ABORT
    162       ENDIF
    163      
    164      
    165       ierr=nf90_inq_varid(nid, string, nvarid)
    166       if (ierr.ne.nf90_noerr) then
    167         write(*,*) 'surfini error, cannot find ',trim(string)
    168         write(*,*) ' in file ',trim(zedatafile),'/surface.nc'
    169         write(*,*)trim(nf90_strerror(ierr))
    170         stop
    171       endif
    172 
    173       ierr=nf90_get_var(nid, nvarid, zdata)
    174 
    175       if (ierr.ne.nf90_noerr) then
    176         write(*,*) 'surfini: error failed loading ',trim(string)
    177         write(*,*)trim(nf90_strerror(ierr))
    178         stop
    179       endif
     192         ENDIF
     193     
     194     
     195         ierr=nf90_inq_varid(nid, string, nvarid)
     196         if (ierr.ne.nf90_noerr) then
     197          write(*,*) 'surfini error, cannot find ',trim(string)
     198          write(*,*) ' in file ',trim(zedatafile),'/surface.nc'
     199          write(*,*)trim(nf90_strerror(ierr))
     200          stop
     201         endif
     202
     203         ierr=nf90_get_var(nid, nvarid, zdata)
     204
     205         if (ierr.ne.nf90_noerr) then
     206          write(*,*) 'surfini: error failed loading ',trim(string)
     207          write(*,*)trim(nf90_strerror(ierr))
     208          stop
     209         endif
    180210 
    181211                     
    182       ierr=nf90_close(nid)
     212         ierr=nf90_close(nid)
    183213 
    184214
    185       nb_ice(:,1) = 1 ! default: there is no ice
    186       latice(:,1) = 1
    187       lonice(:,1) = 1
    188       nb_ice(:,2) = 0
    189       latice(:,2) = 0
    190       lonice(:,2) = 0
    191       !print*,'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
    192 
    193       ! loop over the GCM grid - except for poles (ig=1 and ngrid)
    194       do ig=2,ngrid-1
    195      
    196         ! loop over the surface file grid     
    197         do i=1,imd
    198           do j=1,jmd
    199             zelon = i - 180.
    200             zelat = 90. - j
    201             if ((abs(lati(ig)*180./pi-zelat) .le. 90./real(jjm)) .and.
    202      &        (abs(long(ig)*180./pi-zelon) .le. 180./real(iim))) then
     215         nb_ice(:,1) = 1 ! default: there is no ice
     216         latice(:,1) = 1
     217         lonice(:,1) = 1
     218         nb_ice(:,2) = 0
     219         latice(:,2) = 0
     220         lonice(:,2) = 0
     221         !print*,'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
     222
     223         ! loop over the GCM grid - except for poles (ig=1 and ngrid)
     224         do ig=2,klon_glo-1
     225     
     226          ! loop over the surface file grid     
     227          do i=1,imd
     228           do j=1,jmd
     229             zelon = i - 180.
     230             zelat = 90. - j
     231            if ((abs(lati_glo(ig)*180./pi-zelat).le.90./real(jjm)) .and.
     232     &        (abs(long_glo(ig)*180./pi-zelon).le.180./real(iim))) then
    203233              ! count all points in that GCM grid point
    204234              nb_ice(ig,1) = nb_ice(ig,1) + 1
    … …  
    207237     &           nb_ice(ig,2) = nb_ice(ig,2) + 1
    208238             endif
    209           enddo
    210         enddo 
     239           enddo
     240          enddo 
    211241
    212242        ! projection of nb_ice on GCM lat and lon axes
    213         latice(1+(ig-2)/iim,:) =
     243          latice(1+(ig-2)/iim,:) =
    214244     &     latice(1+(ig-2)/iim,:) + nb_ice(ig,:)
    215         lonice(1+mod(ig-2,iim),:) =
     245          lonice(1+mod(ig-2,iim),:) =
    216246     &     lonice(1+mod(ig-2,iim),:) + nb_ice(ig,:) ! lonice is USELESS ...
    217247
    218       enddo ! of do ig=2,ngrid-1
     248         enddo ! of do ig=2,klon_glo-1
    219249     
    220250
    221251     
    222       ! special case for poles
    223       nb_ice(1,2)   = 1  ! ice prescribed on north pole
    224       latice(1,:)   = nb_ice(1,:)
    225       lonice(1,:)   = nb_ice(1,:)
    226       latice(jjm,:) = nb_ice(ngrid,:)
    227       lonice(iim,:) = nb_ice(ngrid,:)
     252         ! special case for poles
     253         nb_ice(1,2)   = 1  ! ice prescribed on north pole
     254         latice(1,:)   = nb_ice(1,:)
     255         lonice(1,:)   = nb_ice(1,:)
     256         latice(jjm,:) = nb_ice(ngrid,:)
     257         lonice(iim,:) = nb_ice(ngrid,:)
    228258     
    229259     
    … …  
    241271     
    242272   
    243       ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
    244       do i=1,jjm/2
    245       step  = 1. ! threshold to add ice cap
    246       count = 0. ! number of ice GCM caps at this latitude
    247       ! ratiolat is the ratio of area covered by ice within this GCM latitude range
    248       ratiolat  = real(latice(i,2))/real(latice(i,1))
    249       !print*,'i',i,(i-1)*iim+2,i*iim+1
     273         ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
     274         do i=1,jjm/2
     275          step  = 1. ! threshold to add ice cap
     276          count = 0. ! number of ice GCM caps at this latitude
     277          ! ratiolat is the ratio of area covered by ice within this GCM latitude range
     278          ratiolat  = real(latice(i,2))/real(latice(i,1))
     279          !print*,'i',i,(i-1)*iim+2,i*iim+1
    250280     
    251         ! put ice caps while there is not enough ice,
    252         ! as long as the threshold is above 20%
    253         do while ( (count .le. ratiolat*iim ) .and. (step .ge. 0.2))
    254           count = 0.
    255           ! loop over GCM longitudes
    256           do j=1,iim
     281          ! put ice caps while there is not enough ice,
     282          ! as long as the threshold is above 20%
     283          do while ( (count .le. ratiolat*iim ) .and. (step .ge. 0.2))
     284           count = 0.
     285           ! loop over GCM longitudes
     286           do j=1,iim
    257287            ! if the detected ice ratio in the GCM grid point
    258288            ! is more than 'step', then add ice
    259289            if (real(nb_ice((i-1)*iim+1+j,2))
    260290     &        / real(nb_ice((i-1)*iim+1+j,1)) .ge. step) then
    261                   watercaptag((i-1)*iim+1+j) = .true.
     291                  watercaptag_glo((i-1)*iim+1+j) = .true.
    262292                  count = count + 1
    263293            endif
    264           enddo ! of do j=1,iim
     294           enddo ! of do j=1,iim
     295           !print*, 'step',step,count,ratiolat*iim
     296           step = step - 0.01
     297          enddo ! of do while
    265298          !print*, 'step',step,count,ratiolat*iim
    266           step = step - 0.01
    267         enddo ! of do while
    268       !print*, 'step',step,count,ratiolat*iim
    269 
    270       enddo ! of do i=1,jjm/2
     299
     300         enddo ! of do i=1,jjm/2
    271301           
    272302
    273       ELSE IF (icelocationmode .eq. 2) THEN
    274      
    275         print*,'Surfini: predefined ice caps'
    276      
    277         if ((iim .eq. 32) .and. (jjm .eq. 24)) then ! 32x24
     303        ELSE IF (icelocationmode .eq. 2) THEN
     304     
     305         print*,'Surfini: predefined ice caps'
     306     
     307         if ((iim .eq. 32) .and. (jjm .eq. 24)) then ! 32x24
    278308           
    279309          print*,'water ice caps distribution for 32x24 resolution'
    … …  
    288318!---------------------   OUTLIERS  ----------------------------
    289319
    290         else if ((iim .eq. 64) .and. (jjm .eq. 48)) then ! 64x48
     320         else if ((iim .eq. 64) .and. (jjm .eq. 48)) then ! 64x48
    291321
    292322          print*,'water ice caps distribution for 64x48 resolution'
    … …  
    323353
    324354           
    325         else if (ngrid .ne. 1) then
     355         else if (klon_glo .ne. 1) then
    326356       
    327       print*,'No predefined ice location for this resolution :',iim,jjm
    328       print*,'Please change icelocationmode in surfini.F'
    329       print*,'Or add some new definitions ...'
    330       call abort
     357          print*,'No predefined ice location for this resolution :',
     358     &           iim,jjm
     359          print*,'Please change icelocationmode in surfini.F'
     360          print*,'Or add some new definitions ...'
     361          call abort
    331362         
    332         endif
    333 
    334         do ig=1,ngrid
    335           if (longwatercaptag(ig)) watercaptag(ig) = .true.
    336         enddo
    337 
    338 
    339       ELSE IF (icelocationmode .eq. 3) THEN
    340      
    341         print*,'Surfini: ice caps defined by lat and lon values'
    342 
    343          do ig=1,ngrid
     363         endif
     364
     365         do ig=1,klon_glo
     366          if (longwatercaptag(ig)) watercaptag_glo(ig) = .true.
     367         enddo
     368
     369
     370        ELSE IF (icelocationmode .eq. 3) THEN
     371     
     372         print*,'Surfini: ice caps defined by lat and lon values'
     373
     374         do ig=1,klon_glo
    344375         
    345376c-------- Towards olympia planitia water caps -----------
    346377c--------------------------------------------------------
    347378
    348        if ( ( ( lati(ig)*180./pi .ge. 77.  ) .and. ! cap #2
    349      .           ( lati(ig)*180./pi .le. 80.  ) .and.
    350      .           ( long(ig)*180./pi .ge. 110. ) .and.
    351      .           ( long(ig)*180./pi .le. 181. ) )
     379          if ( ( ( lati_glo(ig)*180./pi .ge. 77.  ) .and. ! cap #2
     380     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
     381     .           ( long_glo(ig)*180./pi .ge. 110. ) .and.
     382     .           ( long_glo(ig)*180./pi .le. 181. ) )
    352383     .         .or.
    353384
    354      .         ( ( lati(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
    355      .           ( lati(ig)*180./pi .le. 76.  ) .and.
    356      .           ( long(ig)*180./pi .ge. 150. ) .and.
    357      .           ( long(ig)*180./pi .le. 168. ) )
     385     .         ( ( lati_glo(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
     386     .           ( lati_glo(ig)*180./pi .le. 76.  ) .and.
     387     .           ( long_glo(ig)*180./pi .ge. 150. ) .and.
     388     .           ( long_glo(ig)*180./pi .le. 168. ) )
    358389     .         .or.
    359      .         ( ( lati(ig)*180./pi .ge. 77 ) .and. ! cap #5
    360      .           ( lati(ig)*180./pi .le. 80.  ) .and.
    361      .           ( long(ig)*180./pi .ge. -150.) .and.
    362      .           ( long(ig)*180./pi .le. -110.) ) )
     390     .         ( ( lati_glo(ig)*180./pi .ge. 77 ) .and. ! cap #5
     391     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
     392     .           ( long_glo(ig)*180./pi .ge. -150.) .and.
     393     .           ( long_glo(ig)*180./pi .le. -110.) ) )
    363394     .         then
    364395             
    … …  
    370401               endif !end if alternate = 0
    371402               
    372        endif
     403          endif
    373404
    374405c----------- Opposite olympia planitia water cap --------
    375406c--------------------------------------------------------
    376407
    377         if ( ( ( lati(ig)*180./pi     .ge.  80 ) .and.
    378      .         ( lati(ig)*180./pi     .le.  84 ) )
     408          if ( ( ( lati_glo(ig)*180./pi     .ge.  80 ) .and.
     409     .         ( lati_glo(ig)*180./pi     .le.  84 ) )
    379410     .         .and.
    380      .       ( ( long(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
    381      .         ( long(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
    382 !     .     ( ( ( long(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
    383 !     .         ( long(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
    384 !     .       ( ( long(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
    385 !     .         ( long(ig)*180./pi .le. -70. ) ) ) ) then  !!! 64x48
    386         !   watercaptag(ig)=.true.
    387         endif
     411     .       ( ( long_glo(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
     412     .         ( long_glo(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
     413!     .     ( ( ( long_glo(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
     414!     .         ( long_glo(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
     415!     .       ( ( long_glo(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
     416!     .         ( long_glo(ig)*180./pi .le. -70. ) ) ) ) then  !!! 64x48
     417        !   watercaptag_glo(ig)=.true.
     418          endif
    388419
    389420
    … …  
    391422c--------------------------------------------------------
    392423
    393       if (abs(lati(ig)*180./pi).gt.80)
    394      .       watercaptag(ig)=.true.
     424          if (abs(lati_glo(ig)*180./pi).gt.80)
     425     .          watercaptag_glo(ig)=.true.
    395426           
    396427c--------------------------------------------------------
    397428c--------------------------------------------------------
    398       end do ! of (ngrid)
    399 
    400 
    401        ELSE
     429         end do ! of (klon_glo)
     430
     431
     432        ELSE
    402433     
    403434         print*, 'In surfini.F, icelocationmode is ', icelocationmode
    … …  
    405436         call abort
    406437
    407        ENDIF ! of if (icelocation)
     438        ENDIF ! of if (icelocation)
    408439       
    409440       
    410441        ! print caps locations - useful for plots too
    411         print*,'latitude | longitude | ig'
    412         do ig=1,ngrid
    413           dryness (ig) = icedryness
    414 
    415           if (watercaptag(ig)) then
    416              print*,'ice water cap', lati(ig)*180./pi,
    417      .              long(ig)*180./pi, ig
     442        print*,'surfini: latitude | longitude | ig'
     443        do ig=1,klon_glo
     444          dryness_glo(ig) = icedryness
     445
     446          if (watercaptag_glo(ig)) then
     447             print*,'surfini: ice water cap', lati_glo(ig)*180./pi,
     448     &              long_glo(ig)*180./pi, ig
    418449          endif
    419450        enddo
    420451       
     452       endif !of if (is_master)
     453       
     454       ! Now scatter fields watercaptag and dryness from master to all
     455       ! (is just a plain copy in serial mode)
     456       call scatter(dryness_glo,dryness)
     457       call scatter(watercaptag_glo,watercaptag)
     458       
    421459#endif     
    422 
     460! end of #else of #ifdef MESOSCALE
    423461       ENDIF ! (caps & water)
    424462       
    … …  
    439477              psolaralb(ig,2)=albedodat(ig)
    440478         END DO
    441          PRINT*,'minimum albedo sans water caps',
    442      s     albedodat(ISMIN(ngrid,albedodat,1))
    443          PRINT*,'maximum albedo sans water caps',
    444      s     albedodat(ISMAX(ngrid,albedodat,1))
     479         PRINT*,'surfini: minimum albedo without water caps',
     480     &          minval(albedodat)
     481         PRINT*,'surfini: maximum albedo without water caps',
     482     &          maxval(albedodat)
    445483
    446484c        initial albedo if permanent H2O ice is present
    … …  
    453491            ENDIF
    454492           END DO
    455            PRINT*,'minimum albedo avec water caps',
    456      s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
    457            PRINT*,'maximum albedo avec water caps',
    458      s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
     493           PRINT*,'surfini: minimum albedo with water caps',
     494     &            minval(albedodat)
     495           PRINT*,'surfini: maximum albedo with water caps',
     496     &            maxval(albedodat)
    459497
    460498         ENDIF
    … …  
    465503       DO ig=1,ngrid
    466504         IF (piceco2(ig) .GT. 0.) THEN
    467              IF(ig.GT.ngrid/2+1) THEN
    468                 icap=2
     505             IF(lati(ig).LT. 0.) THEN
     506                icap=2 ! Southern hemisphere
    469507             ELSE
    470                 icap=1
     508                icap=1 ! Northern hemisphere
    471509             ENDIF
    472510             psolaralb(ig,1) = albedice(icap)
    … …  
    514552           endif
    515553          end do
    516           PRINT*,'minimum albedo avec givre et co2',
    517      s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
    518           PRINT*,'maximum albedo avec givre et co2',
    519      s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
     554          PRINT*,'surfini: minimum albedo with frost and co2',
     555     &            minval(albedodat)
     556          PRINT*,'surfini: maximum albedo with frost and co2',
     557     &            maxval(albedodat)
    520558       ELSE
    521559         watercaptag(:) = .false.
    522        END IF
     560       END IF ! OF IF(water)
    523561         
    524562      RETURN

  • trunk/LMDZ.MARS/libf/phymars/tabfi.F

    r1112 r1130  
    4848     &                     iceradius, dtemisice, iceradius
    4949      use yomaer_h, only: tauvis
     50      use iostart, only: get_var
     51      use mod_phys_lmdz_para, only: is_parallel
    5052      implicit none
    5153 
    … …  
    6769c Arguments
    6870c ---------
    69       INTEGER nid,nvarid,tab0
    70       INTEGER*4 day_ini
    71       INTEGER Lmodif
    72       INTEGER lmax
    73       REAL p_rad,p_omeg,p_g,p_mugaz,p_daysec,time,peri_ls
     71      INTEGER,INTENT(IN) :: nid,tab0
     72      INTEGER*4,INTENT(OUT) :: day_ini
     73      INTEGER,INTENT(IN) :: Lmodif
     74      INTEGER,INTENT(OUT) :: lmax
     75      REAL,INTENT(OUT) :: p_rad,p_omeg,p_g,p_mugaz,p_daysec,time
    7476
    7577c Variables
    7678c ---------
     79      INTEGER :: nvarid
     80      REAL :: peri_ls
    7781      INTEGER length
    7882      parameter (length = 100)
    … …  
    8185      INTEGER size
    8286      CHARACTER modif*20
    83 
    84 c Fonctions DRS et autres
    85 c -----------------------
    86       INTEGER setname, cluvdb, getdat
    87 
     87      LOGICAL :: found
     88
     89      write(*,*)"tabfi: nid=",nid," tab0=",tab0," Lmodif=",Lmodif
     90     
    8891c-----------------------------------------------------------------------
    8992c  Initialization of various physical constants to defaut values (nid = 0 case)
    … …  
    150153c Read 'controle' array
    151154c
    152       ierr = NF_INQ_VARID (nid, "controle", nvarid)
    153       IF (ierr .NE. NF_NOERR) THEN
    154          PRINT*, "Tabfi: Could not find <controle> data"
    155          CALL abort
    156       ENDIF
    157 #ifdef NC_DOUBLE
    158       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
    159 #else
    160       ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
    161 #endif
    162       IF (ierr .NE. NF_NOERR) THEN
    163          PRINT*, "Tabfi: Failed reading <controle> array"
    164          CALL abort
    165       ENDIF
    166 
    167        print*,'tabfi: tab_cntrl',tab_cntrl
     155!      ierr = NF_INQ_VARID (nid, "controle", nvarid)
     156!      IF (ierr .NE. NF_NOERR) THEN
     157!         PRINT*, "Tabfi: Could not find <controle> data"
     158!         CALL abort
     159!      ENDIF
     160!#ifdef NC_DOUBLE
     161!      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
     162!#else
     163!      ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
     164!#endif
     165!      IF (ierr .NE. NF_NOERR) THEN
     166!         PRINT*, "Tabfi: Failed reading <controle> array"
     167!         CALL abort
     168!      ENDIF
     169
     170       call get_var("controle",tab_cntrl,found)
     171       if (.not.found) then
     172         write(*,*)"tabfi: Failed reading <controle> array"
     173         call abort
     174       else
     175         write(*,*)'tabfi: tab_cntrl',tab_cntrl
     176       endif
    168177c
    169178c  Initialization of some physical constants
    … …  
    269278c-----------------------------------------------------------------------
    270279c        Modifications...
     280! NB: Modifying controls should only be done by newstart, and in seq mode
     281      if ((Lmodif.eq.1).and.is_parallel) then
     282        write(*,*) "tabfi: Error modifying tab_control should",
     283     &             " only happen in serial mode (eg: by newstart)"
     284        stop
     285      endif
    271286c-----------------------------------------------------------------------
    272287

  • trunk/LMDZ.MARS/libf/phymars/testphys1d.F

    r1112 r1130  
    33! to use  'getin'
    44      USE ioipsl_getincom, only: getin
    5       use infotrac, only: nqtot, tnom
    6       use comsoil_h, only: volcapa, layer, mlayer, inertiedat
     5      use infotrac, only: nqtot, tname
     6      use comsoil_h, only: volcapa, layer, mlayer, inertiedat, nsoilmx
    77      use comgeomfi_h, only: lati, long, area
    88      use comdiurn_h, only: sinlat
    … …  
    1313      use slope_mod, only: theta_sl, psi_sl
    1414      use yomaer_h, only: tauvis
     15      use control_mod, only: day_step
     16      use phyredem, only: physdem0,physdem1
     17      use comgeomphy, only: initcomgeomphy
    1518      IMPLICIT NONE
    1619
    … …  
    3639#include "dimensions.h"
    3740#include "dimphys.h"
     41      integer, parameter :: ngrid = 1 !(2+(jjm-1)*iim - 1/jjm)
    3842!#include "dimradmars.h"
    3943!#include "comgeomfi.h"
    … …  
    4751!#include "comsaison.h"
    4852!#include "yomaer.h"
    49 #include "control.h"
     53!#include "control.h"
    5054#include "comvert.h"
    5155#include "netcdf.inc"
    … …  
    7074      REAL play(nlayermx)   ! Pressure at the middle of the layers (Pa)
    7175      REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
    72       REAL psurf,tsurf     
     76      REAL psurf,tsurf(1)     
    7377      REAL u(nlayermx),v(nlayermx)  ! zonal, meridional wind
    7478      REAL gru,grv   ! prescribed "geostrophic" background wind
    … …  
    7781      REAL,ALLOCATABLE :: qsurf(:) ! tracer surface budget (e.g. kg.m-2)
    7882      REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
    79       REAL co2ice           ! co2ice layer (kg.m-2)
    80       REAL emis             ! surface layer
     83      REAL co2ice(1)        ! co2ice layer (kg.m-2)
     84      REAL emis(1)          ! surface layer
    8185      REAL q2(nlayermx+1)   ! Turbulent Kinetic Energy
    8286      REAL zlay(nlayermx)   ! altitude estimee dans les couches (km)
    … …  
    111115c INITIALISATION
    112116c=======================================================================
     117! initialize "serial/parallel" related stuff
     118!      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
     119      CALL init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))
     120      call initcomgeomphy
    113121
    114122c ------------------------------------------------------
    … …  
    203211        endif
    204212        ! allocate arrays:
    205         allocate(tnom(nq))
     213        allocate(tname(nq))
    206214        allocate(q(nlayermx,nq))
    207215        allocate(qsurf(nq))
    … …  
    212220        ! read tracer names from file traceur.def
    213221        do iq=1,nq
    214           read(90,*,iostat=ierr) tnom(iq)
     222          read(90,*,iostat=ierr) tname(iq)
    215223          if (ierr.ne.0) then
    216224            write(*,*) 'testphys1d: error reading tracer names...'
    … …  
    228236        do iq=1,nq
    229237          txt=""
    230           write(txt,"(a)") tnom(iq)
     238          write(txt,"(a)") tname(iq)
    231239          write(*,*)"  tracer:",trim(txt)
    232240          ! CO2
    … …  
    366374        nq=1
    367375        ! allocate arrays:
    368         allocate(tnom(nq))
     376        allocate(tname(nq))
    369377        allocate(q(nlayermx,nq))
    370378        allocate(qsurf(nq))
    … …  
    374382        do iq=1,nq
    375383          write(str7,'(a1,i2.2)')'q',iq
    376           tnom(iq)=str7
     384          tname(iq)=str7
    377385        enddo
    378386      ! and just to be clean, also initialize tracers to zero for physdem1
    … …  
    549557c  CO2 ice on the surface
    550558c  -------------------
    551       co2ice=0.E+0 ! default value for co2ice
     559      co2ice(1)=0.E+0 ! default value for co2ice
    552560      PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
    553561      call getin("co2ice",co2ice)
    … …  
    558566c  ----------
    559567      emis=emissiv
    560       IF (co2ice.eq.1.E+0) THEN
     568      IF (co2ice(1).eq.1.E+0) THEN
    561569         emis=emisice(1) ! northern hemisphere
    562570         IF(lati(1).LT.0) emis=emisice(2) ! southern hemisphere
    … …  
    615623      DO isoil=1,nsoil
    616624         inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
    617          tsoil(isoil)=tsurf  ! soil temperature
     625         tsoil(isoil)=tsurf(1)  ! soil temperature
    618626      ENDDO
    619627
    … …  
    817825
    818826#include "../dyn3d/disvert.F"
     827#include "../dyn3d/abort_gcm.F"

  • trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90

    r1127 r1130  
    4343
    4444      USE comtherm_h
     45      use planetwide_mod, only: planetwide_maxval
    4546
    4647      IMPLICIT NONE
    … …  
    764765! ===========================================================================
    765766
    766       zlmax=MAXVAL(lmax(:))+2
     767      !zlmax=MAXVAL(lmax(:))+2 ! OK, but in serial mode only; use planet
     768      call planetwide_maxval(lmax,zlmax)
     769      zlmax=zlmax+2
     770     
    767771      if (zlmax .ge. nlayer) then
    768772        print*,'thermals have reached last layer of the model'
    … …  
    790794
    791795! llmax is the index of the heighest thermal in the simulation domain
    792       llmax=1
    793       do ig=1,ngrid
    794          if (lalim(ig)>llmax) llmax=lalim(ig)
    795       enddo
     796      !llmax=1
     797      !do ig=1,ngrid
     798      !   if (lalim(ig)>llmax) llmax=lalim(ig)
     799      !enddo
     800      call planetwide_maxval(lalim,llmax)
    796801
    797802! Integral of a**2/(rho* Delta z), see equation 13 of appendix 4.2 in paper

  • trunk/LMDZ.MARS/libf/phymars/vdif_kc.F

    r1047 r1130  
    3232c
    3333c.......................................................................
    34       REAL dt,g
    35       REAL zlev(ngrid,nlay+1)
    36       REAL zlay(ngrid,nlay)
    37       REAL u(ngrid,nlay)
    38       REAL v(ngrid,nlay)
    39       REAL teta(ngrid,nlay)
    40       REAL cd(ngrid)
    41       REAL q2(ngrid,nlay+1)
    42       REAL km(ngrid,nlay+1)
    43       REAL kn(ngrid,nlay+1)
    44       REAL zq(ngrid,nlay,nq)
     34      REAL,INTENT(IN) :: dt,g
     35      REAL,INTENT(IN) :: zlev(ngrid,nlay+1)
     36      REAL,INTENT(IN) :: zlay(ngrid,nlay)
     37      REAL,INTENT(IN) :: u(ngrid,nlay)
     38      REAL,INTENT(IN) :: v(ngrid,nlay)
     39      REAL,INTENT(IN) :: teta(ngrid,nlay)
     40      REAL,INTENT(IN) :: cd(ngrid)
     41      REAL,INTENT(INOUT) :: q2(ngrid,nlay+1)
     42      REAL,INTENT(OUT) :: km(ngrid,nlay+1)
     43      REAL,INTENT(OUT) :: kn(ngrid,nlay+1)
     44      REAL,INTENT(IN) :: zq(ngrid,nlay,nq)
    4545c.......................................................................
    4646c
    … …  
    247247      nlev=nlay+1
    248248     
     249      ! Other initializations of local variables (to be clean)
     250      long(:,:)=0.
     251      n2(:,:)=0.
     252      sn(:,:)=0.
     253      snq2(:,:)=0.
     254      sm(:,:)=0.
     255      smq2(:,:)=0.
    249256c.......................................................................
    250257c  Special treatment for co2
    … …  
    381388        gninf=.false.
    382389        gnsup=.false.
    383         long(igrid,ilev)=long(igrid,ilev)
    384         long(igrid,ilev)=long(igrid,ilev)
     390        long(igrid,ilev)=long(igrid,ilev) ! not very useful...
     391        long(igrid,ilev)=long(igrid,ilev) ! not very useful...
    385392c
    386393        IF (gn.lt.gnmin) THEN
    … …  
    609616c.......................................................................
    610617c
     618
     619!        call writediagfi(ngrid,'vdif_kc_q2','','',3,q2(:,1:nlay))
     620!        call writediagfi(ngrid,'vdif_kc_km','','',3,km(:,1:nlay))
     621!        call writediagfi(ngrid,'vdif_kc_kn','','',3,kn(:,1:nlay))
     622!        call writediagfi(ngrid,'vdif_kc_unsdz','','',3,unsdz(:,1:nlay))
     623!        call writediagfi(ngrid,'vdif_kc_unsddecz','','',3,
     624!     &                   unsdzdec(:,1:nlay))
     625!        call writediagfi(ngrid,'vdif_kc_q','','',3,q(:,1:nlay))
     626!        call writediagfi(ngrid,'vdif_kc_kmpre','','',3,
     627!     &                   kmpre(:,1:nlay))
     628!        call writediagfi(ngrid,'vdif_kc_long','','',3,long(:,1:nlay))
     629!        call writediagfi(ngrid,'vdif_kc_sn','','',3,sn(:,1:nlay))
     630!        call writediagfi(ngrid,'vdif_kc_sm','','',3,sm(:,1:nlay))
     631
    611632      RETURN
    612633      END

  • trunk/LMDZ.MARS/libf/phymars/vdifc.F

    r1047 r1130  
    5656      REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
    5757      REAL,INTENT(IN) :: ph(ngrid,nlay)
    58       REAL :: pt(ngrid,nlay)
    5958      REAL,INTENT(IN) :: ptsrf(ngrid),pemis(ngrid)
    6059      REAL,INTENT(IN) :: pdufi(ngrid,nlay),pdvfi(ngrid,nlay)
    … …  
    6362      REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
    6463      REAL,INTENT(OUT) :: pdtsrf(ngrid),pdhdif(ngrid,nlay)
    65       REAL,INTENT(IN) ::pcapcal(ngrid)
    66       REAL pq2(ngrid,nlay+1)
     64      REAL,INTENT(IN) :: pcapcal(ngrid)
     65      REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
    6766
    6867c    Argument added for condensation:
    … …  
    8685c   ------
    8786
     87      REAL :: pt(ngrid,nlay)
    8888      REAL ust(ngrid),tst(ngrid)
    8989 
    … …  
    363363! Some usefull diagnostics for the new surface layer parametrization :
    364364
    365 !        call WRITEDIAGFI(ngrid,'pcdv',
     365!        call WRITEDIAGFI(ngrid,'vdifc_zcdv_true',
    366366!     &              'momentum drag','no units',
    367367!     &                         2,zcdv_true)
    368 !        call WRITEDIAGFI(ngrid,'pcdh',
     368!        call WRITEDIAGFI(ngrid,'vdifc_zcdh_true',
    369369!     &              'heat drag','no units',
    370370!     &                         2,zcdh_true)
    371 !        call WRITEDIAGFI(ngrid,'ust',
     371!        call WRITEDIAGFI(ngrid,'vdifc_ust',
    372372!     &              'friction velocity','m/s',2,ust)
    373 !       call WRITEDIAGFI(ngrid,'tst',
     373!       call WRITEDIAGFI(ngrid,'vdifc_tst',
    374374!     &              'friction temperature','K',2,tst)
    375 !        call WRITEDIAGFI(ngrid,'rm-1',
     375!        call WRITEDIAGFI(ngrid,'vdifc_zcdv',
    376376!     &              'aerodyn momentum conductance','m/s',
    377377!     &                         2,zcdv)
    378 !        call WRITEDIAGFI(ngrid,'rh-1',
     378!        call WRITEDIAGFI(ngrid,'vdifc_zcdh',
    379379!     &              'aerodyn heat conductance','m/s',
    380380!     &                         2,zcdh)
    … …  
    384384       IF (.not. calltherm) THEN
    385385
    386        CALL vdif_kc(ptimestep,g,pzlev,pzlay
     386       CALL vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay
    387387     &              ,pu,pv,ph,zcdv_true
    388388     &              ,pq2,zkv,zkh,zq)

  • trunk/LMDZ.MARS/libf/phymars/writediagfi.F

    r1047 r1130  
    4040!=================================================================
    4141      use surfdat_h, only: phisfi
     42      use control_mod, only: ecritphy, day_step, iphysiq
     43      USE mod_phys_lmdz_para, only : is_parallel, is_mpi_root,
     44     &                               is_master, gather
     45      USE mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
    4246      implicit none
    4347
    … …  
    4650!#include "dimphys.h"
    4751#include "paramet.h"
    48 #include "control.h"
     52!#include "control.h"
    4953#include "comvert.h"
    5054#include "comgeom.h"
    … …  
    97101     
    98102#ifndef MESOSCALE
     103
     104#ifdef CPP_PARA
     105! Added to work in parallel mode
     106      real dx3_glop(klon_glo,llm)
     107      real dx3_glo(iim,jjp1,llm) ! to store a global 3D data set
     108      real dx2_glop(klon_glo)
     109      real dx2_glo(iim,jjp1)     ! to store a global 2D (surface) data set
     110      real px2(ngrid)
     111!      real dx1_glo(llm)          ! to store a 1D (column) data set
     112!      real dx0_glo
     113      real phisfi_glo(klon_glo) ! surface geopotential on global physics grid
     114#else
     115      real phisfi_glo(ngrid) ! surface geopotential on global physics grid
     116#endif
     117
    99118!***************************************************************
    100119!Sortie des variables au rythme voulu
    … …  
    162181           stop
    163182         endif
     183         
     184         zitau = -1 ! initialize zitau
     185
     186#ifdef CPP_PARA
     187          ! Gather phisfi() geopotential on physics grid
     188          call Gather(phisfi,phisfi_glo)
     189#else
     190         phisfi_glo(:)=phisfi(:)
     191#endif
     192
     193         !! parallel: we cannot use the usual writediagfi method
     194!!         call iophys_ini
     195         if (is_master) then
     196         ! only the master is required to do this
    164197
    165198         ! Create the NetCDF file
    166          ierr = NF_CREATE(fichnom, IOR(NF_CLOBBER,NF_64BIT_OFFSET), nid)
     199         ierr = NF_CREATE(fichnom, NF_CLOBBER, nid)
    167200         ! Define the 'Time' dimension
    168201         ierr = nf_def_dim(nid,"Time",NF_UNLIMITED,idim)
    … …  
    183216
    184217         ! write "header" of file (longitudes, latitudes, geopotential, ...)
    185          call gr_fi_dyn(1,ngrid,iip1,jjp1,phisfi,phis)
     218         call gr_fi_dyn(1,size(phisfi_glo),iip1,jjp1,phisfi_glo,phis)
    186219         call iniwrite(nid,day_ini,phis)
    187          
    188          zitau = -1 ! initialize zitau
     220
     221         endif ! of if (is_master)
     222
    189223      else
    190          ! Open the NetCDF file
    191          ierr = NF_OPEN(fichnom,NF_WRITE,nid)
     224
     225         if (is_master) then
     226           ! only the master is required to do this
     227
     228           ! Open the NetCDF file
     229           ierr = NF_OPEN(fichnom,NF_WRITE,nid)
     230         endif ! of if (is_master)
     231
    192232      endif ! if (firstnom.eq.'1234567890')
    193233
    … …  
    218258!--------------------------------------------------------
    219259
     260        if (is_master) then
     261           ! only the master is required to do this
    220262        if (nom.eq.firstnom) then
    221263        ! We have identified a "first call" (at given date)
    … …  
    241283        end if ! of if (nom.eq.firstnom)
    242284
     285        endif ! of if (is_master)
     286
    243287!Case of a 3D variable
    244288!---------------------
    245289        if (dim.eq.3) then
    246290
     291#ifdef CPP_PARA
     292          ! Gather field on a "global" (without redundant longitude) array
     293          call Gather(px,dx3_glop)
     294!$OMP MASTER
     295          if (is_mpi_root) then
     296            call Grid1Dto2D_glo(dx3_glop,dx3_glo)
     297            ! copy dx3_glo() to dx3(:) and add redundant longitude
     298            dx3(1:iim,:,:)=dx3_glo(1:iim,:,:)
     299            dx3(iip1,:,:)=dx3(1,:,:)
     300          endif
     301!$OMP END MASTER
     302!$OMP BARRIER
     303#else
    247304!         Passage variable physique -->  variable dynamique
    248305!         recast (copy) variable from physics grid to dynamics grid
    … …  
    260317             ENDDO
    261318           ENDDO
    262 
     319#endif
    263320!         Ecriture du champs
    264321
    265 !         write (*,*) 'In  writediagfi, on sauve:  ' , nom
    266 !         write (*,*) 'In  writediagfi. Estimated date = ' ,date
     322          if (is_master) then
     323           ! only the master writes to output
    267324! name of the variable
    268325           ierr= NF_INQ_VARID(nid,nom,varid)
    … …  
    294351!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
    295352!#else
     353!           write(*,*)"test:  nid=",nid," varid=",varid
     354!           write(*,*)"       corner()=",corner
     355!           write(*,*)"       edges()=",edges
     356!           write(*,*)"       dx3()=",dx3
    296357           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
    297358!#endif
    … …  
    300361              write(*,*) "***** PUT_VAR problem in writediagfi"
    301362              write(*,*) "***** with ",nom
    302               write(*,*) 'ierr=', ierr
     363              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
    303364c             call abort
    304365           endif
    305366
     367          endif !of if (is_master)
     368
    306369!Case of a 2D variable
    307370!---------------------
    308371
    309372        else if (dim.eq.2) then
     373
     374#ifdef CPP_PARA
     375          ! Gather field on a "global" (without redundant longitude) array
     376          px2(:)=px(:,1)
     377          call Gather(px2,dx2_glop)
     378!$OMP MASTER
     379          if (is_mpi_root) then
     380            call Grid1Dto2D_glo(dx2_glop,dx2_glo)
     381            ! copy dx2_glo() to dx2(:) and add redundant longitude
     382            dx2(1:iim,:)=dx2_glo(1:iim,:)
     383            dx2(iip1,:)=dx2(1,:)
     384          endif
     385!$OMP END MASTER
     386!$OMP BARRIER
     387#else
    310388
    311389!         Passage variable physique -->  physique dynamique
    … …  
    323401                dx2(iip1,j)=dx2(1,j)
    324402             ENDDO
    325 
     403#endif
     404
     405          if (is_master) then
     406           ! only the master writes to output
    326407!         write (*,*) 'In  writediagfi, on sauve:  ' , nom
    327408!         write (*,*) 'In  writediagfi. Estimated date = ' ,date
    … …  
    359440              write(*,*) "***** PUT_VAR matter in writediagfi"
    360441              write(*,*) "***** with ",nom
    361               write(*,*) 'ierr=', ierr
     442              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
    362443c             call abort
    363444           endif
    364445
     446          endif !of if (is_master)
     447
    365448!Case of a 1D variable (ie: a column)
    366449!---------------------------------------------------
    367450
    368451       else if (dim.eq.1) then
     452         if (is_parallel) then
     453           write(*,*) "writediagfi error: dim=1 not implemented ",
     454     &                 "in parallel mode"
     455           stop
     456         endif
    369457!         Passage variable physique -->  physique dynamique
    370458!         recast (copy) variable from physics grid to dynamics grid
    … …  
    402490              write(*,*) "***** PUT_VAR problem in writediagfi"
    403491              write(*,*) "***** with ",nom
    404               write(*,*) 'ierr=', ierr
     492              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
    405493c             call abort
    406494           endif
    … …  
    410498
    411499        else if (dim.eq.0) then
     500
    412501           dx0 = px (1,1)
    413502
     503          if (is_master) then
     504           ! only the master writes to output
    414505           ierr= NF_INQ_VARID(nid,nom,varid)
    415506           if (ierr /= NF_NOERR) then
    … …  
    437528              write(*,*) "***** PUT_VAR matter in writediagfi"
    438529              write(*,*) "***** with ",nom
    439               write(*,*) 'ierr=', ierr
     530              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
    440531c             call abort
    441532           endif
    442533
     534          endif !of if (is_master)
     535
    443536        endif ! of if (dim.eq.3) elseif(dim.eq.2)...
    444537
    445538      endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
    446539
    447       ierr= NF_CLOSE(nid)
     540      if (is_master) then
     541        ierr= NF_CLOSE(nid)
     542      endif
    448543
    449544#endif

  • trunk/LMDZ.MARS/libf/phymars/writediagsoil.F90

    r1047 r1130  
    1313
    1414use comsoil_h, only: nsoilmx
     15use control_mod, only: ecritphy, day_step, iphysiq
     16use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
     17use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
    1518
    1619implicit none
    … …  
    1922!#include"dimphys.h"
    2023#include"paramet.h"
    21 #include"control.h"
     24!#include"control.h"
    2225!#include"comsoil.h"
    2326#include"netcdf.inc"
    … …  
    5659integer,dimension(4) :: edges,corners
    5760
     61#ifdef CPP_PARA
     62! Added to work in parallel mode
     63real dx3_glop(klon_glo,nsoilmx)
     64real dx3_glo(iim,jjp1,nsoilmx) ! to store a global 3D data set
     65real dx2_glop(klon_glo)
     66real dx2_glo(iim,jjp1)     ! to store a global 2D (surface) data set
     67real px2(ngrid)
     68#endif
     69
    5870! 1. Initialization step
    5971if (firstname.eq."1234567890") then
    … …  
    7587 
    7688  ! Create output NetCDF file
    77   ierr=NF_CREATE(filename,IOR(NF_CLOBBER,NF_64BIT_OFFSET),nid)
    78   if (ierr.ne.NF_NOERR) then
     89  if (is_master) then
     90   ierr=NF_CREATE(filename,NF_CLOBBER,nid)
     91   if (ierr.ne.NF_NOERR) then
    7992    write(*,*)'writediagsoil: Error, failed creating file '//trim(filename)
    8093    stop
    81   endif
    82  
     94   endif
     95  endif ! of if (is_master)
     96
    8397  ! Define dimensions and axis attributes
    8498  call iniwritesoil(nid,ngrid)
    … …  
    89103else
    90104  ! If not an initialization call, simply open the NetCDF file
    91   ierr=NF_OPEN(filename,NF_WRITE,nid)
     105  if (is_master) then
     106   ierr=NF_OPEN(filename,NF_WRITE,nid)
     107  endif
    92108endif ! of if (firstname.eq."1234567890")
    93109
    … …  
    105121  if (name.eq.firstname) then
    106122    ntime=ntime+1
    107     date=real(zitau+1)/real(day_step)
     123    date=float(zitau+1)/float(day_step)
    108124    ! Note: day_step is known from control.h
    109125   
    110     ! Get NetCDF ID for "time"
    111     ierr=NF_INQ_VARID(nid,"time",varid)
    112     ! Add the current value of date to the "time" array
     126    if (is_master) then
     127     ! Get NetCDF ID for "time"
     128     ierr=NF_INQ_VARID(nid,"time",varid)
     129     ! Add the current value of date to the "time" array
    113130!#ifdef NC_DOUBLE
    114 !    ierr=NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
     131!     ierr=NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
    115132!#else
    116     ierr=NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
     133     ierr=NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
    117134!#endif
    118     if (ierr.ne.NF_NOERR) then
     135     if (ierr.ne.NF_NOERR) then
    119136      write(*,*)"writediagsoil: Failed writing date to time variable"
    120137      stop
    121     endif
     138     endif
     139    endif ! of if (is_master)
    122140  endif ! of if (name.eq.firstname)
    123141
    … …  
    125143if (dimpx.eq.3) then ! Case of a 3D variable
    126144  ! A. Recast data along 'dynamics' grid
     145#ifdef CPP_PARA
     146  ! gather field on a "global" (without redundant longitude) array
     147  call Gather(px,dx3_glop)
     148!$OMP MASTER
     149  if (is_mpi_root) then
     150    call Grid1Dto2D_glo(dx3_glop,dx3_glo)
     151    ! copy dx3_glo() to dx3(:) and add redundant longitude
     152    data3(1:iim,:,:)=dx3_glo(1:iim,:,:)
     153    data3(iip1,:,:)=data3(1,:,:)
     154  endif
     155!$OMP END MASTER
     156!$OMP BARRIER
     157#else
    127158  do l=1,nsoilmx
    128159    ! handle the poles
    … …  
    140171    enddo
    141172  enddo
     173#endif
    142174 
    143175  ! B. Write (append) the variable to the NetCDF file
     176  if (is_master) then
    144177  ! B.1. Get the ID of the variable
    145178  ierr=NF_INQ_VARID(nid,name,varid)
    … …  
    179212               " to file "//trim(filename)//" at time",date
    180213  endif
     214  endif ! of if (is_master)
    181215
    182216elseif (dimpx.eq.2) then ! Case of a 2D variable
     217
    183218  ! A. Recast data along 'dynamics' grid
     219#ifdef CPP_PARA
     220  ! gather field on a "global" (without redundant longitude) array
     221  px2(:)=px(:,1)
     222  call Gather(px2,dx2_glop)
     223!$OMP MASTER
     224  if (is_mpi_root) then
     225    call Grid1Dto2D_glo(dx2_glop,dx2_glo)
     226    ! copy dx3_glo() to dx3(:) and add redundant longitude
     227    data2(1:iim,:)=dx2_glo(1:iim,:)
     228    data2(iip1,:)=data2(1,:)
     229  endif
     230!$OMP END MASTER
     231!$OMP BARRIER
     232#else
    184233  ! handle the poles
    185234  do i=1,iip1
    … …  
    195244    data2(iip1,j)=data2(1,j) ! extra (modulo) longitude
    196245  enddo
     246#endif
    197247
    198248  ! B. Write (append) the variable to the NetCDF file
     249  if (is_master) then
    199250  ! B.1. Get the ID of the variable
    200251  ierr=NF_INQ_VARID(nid,name,varid)
    … …  
    231282               " to file "//trim(filename)//" at time",date
    232283  endif
     284  endif ! of if (is_master)
    233285
    234286elseif (dimpx.eq.0) then ! Case of a 0D variable
     287#ifdef CPP_PARA
     288  write(*,*) "writediagsoil: dimps==0 case not implemented in // mode!!"
     289  stop
     290#endif
    235291  ! A. Copy data value
    236292  data0=px(1,1)
    … …  
    270326
    271327! 4. Close the NetCDF file
    272 ierr=NF_CLOSE(nid)
     328if (is_master) then
     329  ierr=NF_CLOSE(nid)
     330endif
    273331
    274332end subroutine writediagsoil

  • trunk/LMDZ.MARS/libf/phymars/wstats.F90

    r719 r1130  
    11subroutine wstats(ngrid,nom,titre,unite,dim,px)
     2
     3use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather, klon_mpi_begin
     4use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
    25
    36implicit none
    … …  
    58#include "dimensions.h"
    69#include "dimphys.h"
     10#include "comconst.h"
    711#include "statto.h"
    812#include "netcdf.inc"
    … …  
    1115character (len=*),intent(in) :: nom,titre,unite
    1216integer,intent(in) :: dim
    13 real, dimension(ngrid,llm),intent(in) :: px
    1417integer,parameter :: iip1=iim+1
    1518integer,parameter :: jjp1=jjm+1
     19real,intent(in) :: px(ngrid,llm)
    1620real, dimension(iip1,jjp1,llm) :: mean3d,sd3d,dx3
    1721real, dimension(iip1,jjp1) :: mean2d,sd2d,dx2
    1822character (len=50) :: namebis
    1923character (len=50), save :: firstvar
    20 integer :: ierr,ierr2,varid,nbdim,nid
     24integer :: ierr,varid,nbdim,nid
    2125integer :: meanid,sdid
    22 integer, dimension(4)  :: id,start,size
     26integer, dimension(4)  :: id,start,sizes
    2327logical, save :: firstcall=.TRUE.
    2428integer :: l,i,j,ig0
    25 integer,save :: index
    26 ! Added to use stats.def to select output variable
    27 logical,save :: stats_def
    28 logical :: getout
    29 integer,save :: n_nom_def
    30 integer :: n
    31 integer,parameter :: n_nom_def_max=99
    32 character(len=20),save :: nom_def(n_nom_def_max)
    33 logical, save :: notfoundfirst=.TRUE.
     29integer,save :: indx
    3430
    3531integer, save :: step=0
    36      
    37      
    38 
     32
     33! Added to work in parallel mode
     34#ifdef CPP_PARA
     35real px3_glop(klon_glo,llm) ! to store a 3D data set on global physics grid
     36real px3_glo(iim,jjp1,llm) ! to store a global 3D data set on global lonxlat grid
     37real px2_glop(klon_glo) ! to store a 2D data set on global physics grid
     38real px2_glo(iim,jjp1) ! to store a 2D data set on global lonxlat grid
     39real px2(ngrid)
     40real px3(ngrid,llm)
     41#else
     42! When not running in parallel mode:
     43real px3_glop(ngrid,llm) ! to store a 3D data set on global physics grid
     44real px3_glo(iim,jjp1,llm) ! to store a global 3D data set on global lonxlat grid
     45real px2_glop(ngrid) ! to store a 2D data set on global physics grid
     46real px2_glo(iim,jjp1) ! to store a 2D data set on global lonxlat grid
     47#endif
     48
     49! 1. Initialization (creation of stats.nc file)
    3950if (firstcall) then
    4051   firstcall=.false.
     52   firstvar=trim((nom))
    4153   call inistats(ierr)
    42 ! at very first call, check if there is a "stats.def" to use and read it
    43    open(99,file="stats.def",status='old',form='formatted',iostat=ierr2)
    44    if(ierr2.eq.0) then
    45      stats_def=.true.
    46      write(*,*) "******************"
    47      write(*,*) "Reading stats.def"
    48      write(*,*) "******************"
    49      do n=1,n_nom_def_max
    50           read(99,fmt='(a)',end=88) nom_def(n)
    51           write(*,*) 'Output in stats: ', trim(nom_def(n))
    52      end do
    53 88   continue
    54      if (n.ge.n_nom_def_max) then
    55         write(*,*)"n_nom_def_max too small in wstats.F90:",n
    56         stop
    57      end if
    58      n_nom_def=n-1
    59      close(99)
    60    else
    61      firstvar=trim((nom))
    62      stats_def=.false.
    63      notfoundfirst=.false.
    64    endif
    65 endif
    66 
    67 ! find firstvar that is in stats.def
    68 if (notfoundfirst) then
    69    do n=1,n_nom_def_max
    70      if(trim((nom_def(n))).eq.trim((nom))) then
    71         firstvar=trim((nom))
    72         notfoundfirst=.false.
    73      endif
    74    end do
    75 endif
    76 
    77 ! get out of wstats if there is stats.def AND variable not listed
    78 if (stats_def) then
    79    getout=.true.
    80    do n=1,n_nom_def
    81        if(trim(nom_def(n)).eq.nom) getout=.false.
    82    end do
    83    if (getout) return
    84 end if
    85 
    86 if (firstvar==nom) then ! If we're back to the first variable
     54endif
     55
     56if (firstvar==nom) then ! If we're back to the first variable, increment time counter
    8757      step = step + 1
    8858endif
    8959
    9060if (mod(step,istats).ne.0) then
     61  ! if its not time to write to file, exit
    9162   RETURN
    9263endif
    9364
     65! collect fields on a global physics grid
     66#ifdef CPP_PARA
     67 if (dim.eq.3) then
     68  px3(1:ngrid,1:llm)=px(1:ngrid,1:llm)
     69  ! Gather fieds on a "global" (without redundant longitude) array
     70  call Gather(px3,px3_glop)
     71!$OMP MASTER
     72  if (is_mpi_root) then
     73    call Grid1Dto2D_glo(px3_glop,px3_glo)
     74    ! copy dx3_glo() to dx3(:) and add redundant longitude
     75    dx3(1:iim,:,:)=px3_glo(1:iim,:,:)
     76    dx3(iip1,:,:)=dx3(1,:,:)
     77  endif
     78!$OMP END MASTER
     79!$OMP BARRIER
     80 else ! dim.eq.2
     81  ! Gather fieds on a "global" (without redundant longitude) array
     82  px2(:)=px(:,1)
     83  call Gather(px2,px2_glop)
     84!$OMP MASTER
     85          if (is_mpi_root) then
     86            call Grid1Dto2D_glo(px2_glop,px2_glo)
     87            ! copy px2_glo() to dx2(:) and add redundant longitude
     88            dx2(1:iim,:)=px2_glo(1:iim,:)
     89            dx2(iip1,:)=dx2(1,:)
     90          endif
     91!$OMP END MASTER
     92!$OMP BARRIER
     93 endif
     94#else
     95  if (dim.eq.3) then
     96    px3_glop(:,1:llm)=px(:,1:llm)
     97!  Passage variable physique -->  variable dynamique
     98    DO l=1,llm
     99      DO i=1,iim
     100         px3_glo(i,1,l)=px(1,l)
     101         px3_glo(i,jjp1,l)=px(ngrid,l)
     102      ENDDO
     103      DO j=2,jjm
     104         ig0= 1+(j-2)*iim
     105         DO i=1,iim
     106            px3_glo(i,j,l)=px(ig0+i,l)
     107         ENDDO
     108      ENDDO
     109    ENDDO
     110  else ! dim.eq.2
     111    px2_glop(:)=px(:,1)
     112!    Passage variable physique -->  physique dynamique
     113   DO i=1,iim
     114     px2_glo(i,1)=px(1,1)
     115     px2_glo(i,jjp1)=px(ngrid,1)
     116   ENDDO
     117   DO j=2,jjm
     118     ig0= 1+(j-2)*iim
     119     DO i=1,iim
     120        px2_glo(i,j)=px(ig0+i,1)
     121     ENDDO
     122   ENDDO
     123  endif
     124#endif
     125
     126! 2. Write field to file
     127
     128if (is_master) then
     129! only master needs do this
     130
    94131ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
    95132
    96133namebis=trim(nom)
     134
     135! test: check if that variable already exists in file
    97136ierr= NF_INQ_VARID(nid,namebis,meanid)
    98137
    99138if (ierr.ne.NF_NOERR) then
    100 
     139  ! variable not in file, create/define it
    101140   if (firstvar==nom) then
    102       index=1
    103       count=0
     141      indx=1
     142      count(:)=0
    104143   endif
    105144
    … …  
    122161   write (*,*) "STATS: creation de ",nom
    123162   namebis=trim(nom)
    124    call def_var_stats(nid,namebis,titre,unite,nbdim,id,meanid,ierr)
    125    call inivar(nid,meanid,ngrid,dim,index,px,ierr)
     163   call def_var(nid,namebis,titre,unite,nbdim,id,meanid,ierr)
     164   if (dim.eq.3) then
     165     call inivar(nid,meanid,size(px3_glop,1),dim,indx,px3_glop,ierr)
     166   else ! dim.eq.2
     167     call inivar(nid,meanid,size(px2_glop,1),dim,indx,px2_glop,ierr)
     168   endif
    126169   namebis=trim(nom)//"_sd"
    127    call def_var_stats(nid,namebis,trim(titre)//" total standard deviation over the season",unite,nbdim,id,sdid,ierr)
    128    call inivar(nid,sdid,ngrid,dim,index,px,ierr)
     170   call def_var(nid,namebis,trim(titre)//" total standard deviation over the season",unite,nbdim,id,sdid,ierr)
     171   if (dim.eq.3) then
     172     call inivar(nid,sdid,size(px3_glop,1),dim,indx,px3_glop,ierr)
     173   else ! dim.eq.2
     174     call inivar(nid,sdid,size(px2_glop,1),dim,indx,px2_glop,ierr)
     175   endif
    129176
    130177   ierr= NF_CLOSE(nid)
    … …  
    132179
    133180else
     181   ! variable found in file
    134182   namebis=trim(nom)//"_sd"
    135183   ierr= NF_INQ_VARID(nid,namebis,sdid)
    … …  
    138186
    139187if (firstvar==nom) then
    140    count(index)=count(int(index))+1
    141    index=index+1
    142    if (index>istime) then
    143       index=1
    144    endif
    145 endif
    146 
    147 if (count(index)==0) then
     188   count(indx)=count(int(indx))+1
     189   indx=indx+1
     190   if (indx>istime) then
     191      indx=1
     192   endif
     193endif
     194
     195if (count(indx)==0) then
     196   ! very first time we write the variable to file
    148197   if (dim.eq.3) then
    149       start=(/1,1,1,index/)
    150       size=(/iip1,jjp1,llm,1/)
    151       mean3d=0
    152       sd3d=0
     198      start=(/1,1,1,indx/)
     199      sizes=(/iip1,jjp1,llm,1/)
     200      mean3d(:,:,:)=0
     201      sd3d(:,:,:)=0
    153202   else if (dim.eq.2) then
    154       start=(/1,1,index,0/)
    155       size=(/iip1,jjp1,1,0/)
    156       mean2d=0
    157       sd2d=0
     203      start=(/1,1,indx,0/)
     204      sizes=(/iip1,jjp1,1,0/)
     205      mean2d(:,:)=0
     206      sd2d(:,:)=0
    158207   endif
    159208else
     209   ! load values from file
    160210   if (dim.eq.3) then
    161       start=(/1,1,1,index/)
    162       size=(/iip1,jjp1,llm,1/)
    163 #ifdef NC_DOUBLE
    164       ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,size,mean3d)
    165       ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,size,sd3d)
    166 #else
    167       ierr = NF_GET_VARA_REAL(nid,meanid,start,size,mean3d)
    168       ierr = NF_GET_VARA_REAL(nid,sdid,start,size,sd3d)
     211      start=(/1,1,1,indx/)
     212      sizes=(/iip1,jjp1,llm,1/)
     213#ifdef NC_DOUBLE
     214      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
     215      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
     216#else
     217      ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean3d)
     218      ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd3d)
    169219#endif
    170220      if (ierr.ne.NF_NOERR) then
    … …  
    174224
    175225   else if (dim.eq.2) then
    176       start=(/1,1,index,0/)
    177       size=(/iip1,jjp1,1,0/)
    178 #ifdef NC_DOUBLE
    179       ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,size,mean2d)
    180       ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,size,sd2d)
    181 #else
    182       ierr = NF_GET_VARA_REAL(nid,meanid,start,size,mean2d)
    183       ierr = NF_GET_VARA_REAL(nid,sdid,start,size,sd2d)
     226      start=(/1,1,indx,0/)
     227      sizes=(/iip1,jjp1,1,0/)
     228#ifdef NC_DOUBLE
     229      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
     230      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
     231#else
     232      ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean2d)
     233      ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd2d)
    184234#endif
    185235      if (ierr.ne.NF_NOERR) then
    … …  
    188238      endif
    189239   endif
    190 endif
     240endif ! of if (count(indx)==0)
     241
     242
     243! 2.1. Build dx* (data on lon-lat grid, with redundant longitude)
    191244
    192245if (dim.eq.3) then
     246  dx3(1:iim,:,:)=px3_glo(:,:,:)
     247  dx3(iip1,:,:)=dx3(1,:,:)
     248else ! dim.eq.2
     249  dx2(1:iim,:)=px2_glo(:,:)
     250  dx2(iip1,:)=dx2(1,:)
     251endif
     252
     253
     254! 2.2. Add current values to previously stored sums
     255
     256if (dim.eq.3) then
     257
     258   mean3d(:,:,:)=mean3d(:,:,:)+dx3(:,:,:)
     259   sd3d(:,:,:)=sd3d(:,:,:)+dx3(:,:,:)**2
     260
     261#ifdef NC_DOUBLE
     262   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
     263   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
     264#else
     265   ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean3d)
     266   ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd3d)
     267#endif
     268
     269else if (dim.eq.2) then
     270
     271   mean2d(:,:)= mean2d(:,:)+dx2(:,:)
     272   sd2d(:,:)=sd2d(:,:)+dx2(:,:)**2
     273
     274#ifdef NC_DOUBLE
     275   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
     276   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
     277#else
     278   ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean2d)
     279   ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd2d)
     280#endif
     281
     282endif ! of if (dim.eq.3) elseif (dim.eq.2)
     283
     284  ierr= NF_CLOSE(nid)
     285endif ! of if (is_master)
     286
     287end subroutine wstats
     288
     289!======================================================
     290subroutine inivar(nid,varid,ngrid,dim,indx,px,ierr)
     291
     292implicit none
     293
     294include "dimensions.h"
     295include "dimphys.h"
     296include "netcdf.inc"
     297
     298integer, intent(in) :: nid,varid,dim,indx,ngrid
     299real, dimension(ngrid,llm), intent(in) :: px
     300integer, intent(out) :: ierr
     301
     302integer,parameter :: iip1=iim+1
     303integer,parameter :: jjp1=jjm+1
     304
     305integer :: l,i,j,ig0
     306integer, dimension(4) :: start,sizes
     307real, dimension(iip1,jjp1,llm) :: dx3
     308real, dimension(iip1,jjp1) :: dx2
     309
     310if (dim.eq.3) then
     311
     312   start=(/1,1,1,indx/)
     313   sizes=(/iip1,jjp1,llm,1/)
    193314
    194315!  Passage variable physique -->  variable dynamique
    … …  
    208329   ENDDO
    209330
    210    mean3d(:,:,:)= mean3d(:,:,:)+dx3(:,:,:)
    211    sd3d(:,:,:)= sd3d(:,:,:)+dx3(:,:,:)**2
    212 
    213 #ifdef NC_DOUBLE
    214    ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,size,mean3d)
    215    ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,size,sd3d)
    216 #else
    217    ierr = NF_PUT_VARA_REAL(nid,meanid,start,size,mean3d)
    218    ierr = NF_PUT_VARA_REAL(nid,sdid,start,size,sd3d)
    219 #endif
    220    if (ierr.ne.NF_NOERR) then
    221      write (*,*) NF_STRERROR(ierr)
    222      stop ""
    223    endif
     331#ifdef NC_DOUBLE
     332   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx3)
     333#else
     334   ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx3)
     335#endif
    224336
    225337else if (dim.eq.2) then
     338
     339      start=(/1,1,indx,0/)
     340      sizes=(/iip1,jjp1,1,0/)
    226341
    227342!    Passage variable physique -->  physique dynamique
    … …  
    239354  ENDDO
    240355
    241    mean2d(:,:)= mean2d(:,:)+dx2(:,:)
    242    sd2d(:,:)= sd2d(:,:)+dx2(:,:)**2
    243 
    244 #ifdef NC_DOUBLE
    245    ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,size,mean2d)
    246    ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,size,sd2d)
    247 #else
    248    ierr = NF_PUT_VARA_REAL(nid,meanid,start,size,mean2d)
    249    ierr = NF_PUT_VARA_REAL(nid,sdid,start,size,sd2d)
    250 #endif
    251    if (ierr.ne.NF_NOERR) then
    252      write (*,*) NF_STRERROR(ierr)
    253      stop ""
    254    endif
    255 
    256 endif
    257 
    258 ierr= NF_CLOSE(nid)
    259 
    260 end subroutine wstats
    261 
    262 !======================================================
    263 subroutine inivar(nid,varid,ngrid,dim,index,px,ierr)
    264 
    265 implicit none
    266 
    267 include "dimensions.h"
    268 include "dimphys.h"
    269 include "netcdf.inc"
    270 
    271 integer, intent(in) :: nid,varid,dim,index,ngrid
    272 real, dimension(ngrid,llm), intent(in) :: px
    273 integer, intent(out) :: ierr
    274 
    275 integer,parameter :: iip1=iim+1
    276 integer,parameter :: jjp1=jjm+1
    277 
    278 integer :: l,i,j,ig0
    279 integer, dimension(4) :: start,size
    280 real, dimension(iip1,jjp1,llm) :: dx3
    281 real, dimension(iip1,jjp1) :: dx2
    282 
    283 if (dim.eq.3) then
    284 
    285    start=(/1,1,1,index/)
    286    size=(/iip1,jjp1,llm,1/)
    287 
    288 !  Passage variable physique -->  variable dynamique
    289 
    290    DO l=1,llm
    291       DO i=1,iip1
    292          dx3(i,1,l)=px(1,l)
    293          dx3(i,jjp1,l)=px(ngrid,l)
    294       ENDDO
    295       DO j=2,jjm
    296          ig0= 1+(j-2)*iim
    297          DO i=1,iim
    298             dx3(i,j,l)=px(ig0+i,l)
    299          ENDDO
    300          dx3(iip1,j,l)=dx3(1,j,l)
    301       ENDDO
    302    ENDDO
    303 
    304 #ifdef NC_DOUBLE
    305    ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,dx3)
    306 #else
    307    ierr = NF_PUT_VARA_REAL(nid,varid,start,size,dx3)
    308 #endif
    309 
    310 else if (dim.eq.2) then
    311 
    312       start=(/1,1,index,0/)
    313       size=(/iip1,jjp1,1,0/)
    314 
    315 !    Passage variable physique -->  physique dynamique
    316 
    317   DO i=1,iip1
    318      dx2(i,1)=px(1,1)
    319      dx2(i,jjp1)=px(ngrid,1)
    320   ENDDO
    321   DO j=2,jjm
    322      ig0= 1+(j-2)*iim
    323      DO i=1,iim
    324         dx2(i,j)=px(ig0+i,1)
    325      ENDDO
    326      dx2(iip1,j)=dx2(1,j)
    327   ENDDO
    328 
    329 #ifdef NC_DOUBLE
    330    ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,dx2)
    331 #else
    332    ierr = NF_PUT_VARA_REAL(nid,varid,start,size,dx2)
     356#ifdef NC_DOUBLE
     357   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx2)
     358#else
     359   ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx2)
    333360#endif
    334361

  • trunk/LMDZ.MARS/libf/phymars/xvik.F

    r1087 r1130  
    1919#include "logic.h"
    2020#include "temps.h"
    21 #include "control.h"
     21!#include "control.h"
    2222#include "ener.h"
    2323#include "description.h"
Note: See TracChangeset for help on using the changeset viewer.