Changeset 1188 for LMDZ4


Ignore:
Timestamp:
Jun 18, 2009, 5:50:11 PM (15 years ago)
Author:
lguez
Message:

Translated calls using NetCDF 2.4 interface to calls using NetCDF 3.6
Fortran 90 interface.
Corrected a few comments in "output.def".

Location:
LMDZ4/branches/LMDZ4-dev
Files:
4 edited
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ4/branches/LMDZ4-dev/libf/dyn3d/guide_mod.F90

    r1187 r1188  
    1111  USE getparam
    1212  USE Write_Field
     13  use netcdf, only: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close
    1314
    1415  IMPLICIT NONE
    … …  
    140141    ncidpl=-99
    141142    if (guide_modele) then
    142        if (ncidpl.eq.-99) ncidpl=NCOPN('apbp.nc',NCNOWRIT,rcod)
     143       if (ncidpl.eq.-99) rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
    143144    else
    144145         if (guide_u) then
    145            if (ncidpl.eq.-99) ncidpl=NCOPN('u.nc',NCNOWRIT,rcod)
     146           if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
    146147         elseif (guide_v) then
    147            if (ncidpl.eq.-99) ncidpl=NCOPN('v.nc',NCNOWRIT,rcod)
     148           if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
    148149         elseif (guide_T) then
    149            if (ncidpl.eq.-99) ncidpl=NCOPN('T.nc',NCNOWRIT,rcod)
     150           if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
    150151         elseif (guide_Q) then
    151            if (ncidpl.eq.-99) ncidpl=NCOPN('hur.nc',NCNOWRIT,rcod)
     152           if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
    152153         endif
    153154    endif
    … …  
    160161    error=NF_INQ_DIMLEN(ncidpl,rid,nlevnc)
    161162    print *,'Guide: nombre niveaux vert. nlevnc', nlevnc
    162     CALL NCCLOS(ncidpl,rcod)
     163    rcod = nf90_close(ncidpl)
    163164
    164165! ---------------------------------------------
    … …  
    992993         if (guide_modele) then
    993994             print *,'Lecture du guidage sur niveaux mod�le'
    994              ncidpl=NCOPN('apbp.nc',NCNOWRIT,rcode)
    995              varidap=NCVID(ncidpl,'AP',rcode)
    996              varidbp=NCVID(ncidpl,'BP',rcode)
     995             rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
     996             rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
     997             rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
    997998             print*,'ncidpl,varidap',ncidpl,varidap
    998999         endif
    9991000! Vent zonal
    10001001         if (guide_u) then
    1001              ncidu=NCOPN('u.nc',NCNOWRIT,rcode)
    1002              varidu=NCVID(ncidu,'UWND',rcode)
     1002             rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
     1003             rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
    10031004             print*,'ncidu,varidu',ncidu,varidu
    10041005             if (ncidpl.eq.-99) ncidpl=ncidu
    … …  
    10061007! Vent meridien
    10071008         if (guide_v) then
    1008              ncidv=NCOPN('v.nc',NCNOWRIT,rcode)
    1009              varidv=NCVID(ncidv,'VWND',rcode)
     1009             rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
     1010             rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
    10101011             print*,'ncidv,varidv',ncidv,varidv
    10111012             if (ncidpl.eq.-99) ncidpl=ncidv
    … …  
    10131014! Temperature
    10141015         if (guide_T) then
    1015              ncidt=NCOPN('T.nc',NCNOWRIT,rcode)
    1016              varidt=NCVID(ncidt,'AIR',rcode)
     1016             rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
     1017             rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
    10171018             print*,'ncidT,varidT',ncidt,varidt
    10181019             if (ncidpl.eq.-99) ncidpl=ncidt
    … …  
    10201021! Humidite
    10211022         if (guide_Q) then
    1022              ncidQ=NCOPN('hur.nc',NCNOWRIT,rcode)
    1023              varidQ=NCVID(ncidQ,'RH',rcode)
     1023             rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
     1024             rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
    10241025             print*,'ncidQ,varidQ',ncidQ,varidQ
    10251026             if (ncidpl.eq.-99) ncidpl=ncidQ
    … …  
    10271028! Pression de surface
    10281029         if ((guide_P).OR.(guide_modele)) then
    1029              ncidps=NCOPN('ps.nc',NCNOWRIT,rcode)
    1030              varidps=NCVID(ncidps,'SP',rcode)
     1030             rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
     1031             rcode = nf90_inq_varid(ncidps, 'SP', varidps)
    10311032             print*,'ncidps,varidps',ncidps,varidps
    10321033         endif
    10331034! Coordonnee verticale
    10341035         if (.not.guide_modele) then
    1035               rcode=NF_INQ_VARID(ncidpl,'LEVEL',varidpl)
    1036               IF (rcode.NE.0) rcode=NF_INQ_VARID(ncidpl,'PRESSURE',varidpl)
     1036              rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
     1037              IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
    10371038              print*,'ncidpl,varidpl',ncidpl,varidpl
    10381039         endif
    … …  
    11751176         if (guide_modele) then
    11761177             print *,'Lecture du guidage sur niveaux mod�le'
    1177              ncidpl=NCOPN('apbp.nc',NCNOWRIT,rcode)
    1178              varidap=NCVID(ncidpl,'AP',rcode)
    1179              varidbp=NCVID(ncidpl,'BP',rcode)
     1178             rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
     1179             rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
     1180             rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
    11801181             print*,'ncidpl,varidap',ncidpl,varidap
    11811182         endif
    11821183! Vent zonal
    11831184         if (guide_u) then
    1184              ncidu=NCOPN('u.nc',NCNOWRIT,rcode)
    1185              varidu=NCVID(ncidu,'UWND',rcode)
     1185             rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
     1186             rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
    11861187             print*,'ncidu,varidu',ncidu,varidu
    11871188             if (ncidpl.eq.-99) ncidpl=ncidu
    … …  
    11891190! Vent meridien
    11901191         if (guide_v) then
    1191              ncidv=NCOPN('v.nc',NCNOWRIT,rcode)
    1192              varidv=NCVID(ncidv,'VWND',rcode)
     1192             rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
     1193             rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
    11931194             print*,'ncidv,varidv',ncidv,varidv
    11941195             if (ncidpl.eq.-99) ncidpl=ncidv
    … …  
    11961197! Temperature
    11971198         if (guide_T) then
    1198              ncidt=NCOPN('T.nc',NCNOWRIT,rcode)
    1199              varidt=NCVID(ncidt,'AIR',rcode)
     1199             rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
     1200             rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
    12001201             print*,'ncidT,varidT',ncidt,varidt
    12011202             if (ncidpl.eq.-99) ncidpl=ncidt
    … …  
    12031204! Humidite
    12041205         if (guide_Q) then
    1205              ncidQ=NCOPN('hur.nc',NCNOWRIT,rcode)
    1206              varidQ=NCVID(ncidQ,'RH',rcode)
     1206             rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
     1207             rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
    12071208             print*,'ncidQ,varidQ',ncidQ,varidQ
    12081209             if (ncidpl.eq.-99) ncidpl=ncidQ
    … …  
    12101211! Pression de surface
    12111212         if ((guide_P).OR.(guide_modele)) then
    1212              ncidps=NCOPN('ps.nc',NCNOWRIT,rcode)
    1213              varidps=NCVID(ncidps,'SP',rcode)
     1213             rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
     1214             rcode = nf90_inq_varid(ncidps, 'SP', varidps)
    12141215             print*,'ncidps,varidps',ncidps,varidps
    12151216         endif
    12161217! Coordonnee verticale
    12171218         if (.not.guide_modele) then
    1218               rcode=NF_INQ_VARID(ncidpl,'LEVEL',varidpl)
    1219               IF (rcode.NE.0) rcode=NF_INQ_VARID(ncidpl,'PRESSURE',varidpl)
     1219              rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
     1220              IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
    12201221              print*,'ncidpl,varidpl',ncidpl,varidpl
    12211222         endif

  • LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/guide_p_mod.F90

    r1187 r1188  
    1212  USE getparam
    1313  USE Write_Field_p
     14  use netcdf, only: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close
    1415
    1516  IMPLICIT NONE
    … …  
    145146    ncidpl=-99
    146147    if (guide_modele) then
    147        if (ncidpl.eq.-99) ncidpl=NCOPN('apbp.nc',NCNOWRIT,rcod)
     148       if (ncidpl.eq.-99) rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
    148149    else
    149150         if (guide_u) then
    150            if (ncidpl.eq.-99) ncidpl=NCOPN('u.nc',NCNOWRIT,rcod)
     151           if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
    151152         elseif (guide_v) then
    152            if (ncidpl.eq.-99) ncidpl=NCOPN('v.nc',NCNOWRIT,rcod)
     153           if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
    153154         elseif (guide_T) then
    154            if (ncidpl.eq.-99) ncidpl=NCOPN('T.nc',NCNOWRIT,rcod)
     155           if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
    155156         elseif (guide_Q) then
    156            if (ncidpl.eq.-99) ncidpl=NCOPN('hur.nc',NCNOWRIT,rcod)
     157           if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
    157158         endif
    158159    endif
    … …  
    165166    error=NF_INQ_DIMLEN(ncidpl,rid,nlevnc)
    166167    print *,'Guide: nombre niveaux vert. nlevnc', nlevnc
    167     CALL NCCLOS(ncidpl,rcod)
     168    rcod = nf90_close(ncidpl)
    168169
    169170! ---------------------------------------------
    … …  
    10711072         if (guide_modele) then
    10721073             print *,'Lecture du guidage sur niveaux mod�le'
    1073              ncidpl=NCOPN('apbp.nc',NCNOWRIT,rcode)
    1074              varidap=NCVID(ncidpl,'AP',rcode)
    1075              varidbp=NCVID(ncidpl,'BP',rcode)
     1074             rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
     1075             rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
     1076             rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
    10761077             print*,'ncidpl,varidap',ncidpl,varidap
    10771078         endif
    10781079! Vent zonal
    10791080         if (guide_u) then
    1080              ncidu=NCOPN('u.nc',NCNOWRIT,rcode)
    1081              varidu=NCVID(ncidu,'UWND',rcode)
     1081             rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
     1082             rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
    10821083             print*,'ncidu,varidu',ncidu,varidu
    10831084             if (ncidpl.eq.-99) ncidpl=ncidu
    … …  
    10851086! Vent meridien
    10861087         if (guide_v) then
    1087              ncidv=NCOPN('v.nc',NCNOWRIT,rcode)
    1088              varidv=NCVID(ncidv,'VWND',rcode)
     1088             rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
     1089             rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
    10891090             print*,'ncidv,varidv',ncidv,varidv
    10901091             if (ncidpl.eq.-99) ncidpl=ncidv
    … …  
    10921093! Temperature
    10931094         if (guide_T) then
    1094              ncidt=NCOPN('T.nc',NCNOWRIT,rcode)
    1095              varidt=NCVID(ncidt,'AIR',rcode)
     1095             rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
     1096             rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
    10961097             print*,'ncidT,varidT',ncidt,varidt
    10971098             if (ncidpl.eq.-99) ncidpl=ncidt
    … …  
    10991100! Humidite
    11001101         if (guide_Q) then
    1101              ncidQ=NCOPN('hur.nc',NCNOWRIT,rcode)
    1102              varidQ=NCVID(ncidQ,'RH',rcode)
     1102             rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
     1103             rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
    11031104             print*,'ncidQ,varidQ',ncidQ,varidQ
    11041105             if (ncidpl.eq.-99) ncidpl=ncidQ
    … …  
    11061107! Pression de surface
    11071108         if ((guide_P).OR.(guide_modele)) then
    1108              ncidps=NCOPN('ps.nc',NCNOWRIT,rcode)
    1109              varidps=NCVID(ncidps,'SP',rcode)
     1109             rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
     1110             rcode = nf90_inq_varid(ncidps, 'SP', varidps)
    11101111             print*,'ncidps,varidps',ncidps,varidps
    11111112         endif
    11121113! Coordonnee verticale
    11131114         if (.not.guide_modele) then
    1114               rcode=NF_INQ_VARID(ncidpl,'LEVEL',varidpl)
    1115               IF (rcode.NE.0) rcode=NF_INQ_VARID(ncidpl,'PRESSURE',varidpl)
     1115              rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
     1116              IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
    11161117              print*,'ncidpl,varidpl',ncidpl,varidpl
    11171118         endif
    … …  
    12551256         if (guide_modele) then
    12561257             print *,'Lecture du guidage sur niveaux mod�le'
    1257              ncidpl=NCOPN('apbp.nc',NCNOWRIT,rcode)
    1258              varidap=NCVID(ncidpl,'AP',rcode)
    1259              varidbp=NCVID(ncidpl,'BP',rcode)
     1258             rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
     1259             rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
     1260             rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
    12601261             print*,'ncidpl,varidap',ncidpl,varidap
    12611262         endif
    12621263! Vent zonal
    12631264         if (guide_u) then
    1264              ncidu=NCOPN('u.nc',NCNOWRIT,rcode)
    1265              varidu=NCVID(ncidu,'UWND',rcode)
     1265             rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
     1266             rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
    12661267             print*,'ncidu,varidu',ncidu,varidu
    12671268             if (ncidpl.eq.-99) ncidpl=ncidu
    … …  
    12691270! Vent meridien
    12701271         if (guide_v) then
    1271              ncidv=NCOPN('v.nc',NCNOWRIT,rcode)
    1272              varidv=NCVID(ncidv,'VWND',rcode)
     1272             rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
     1273             rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
    12731274             print*,'ncidv,varidv',ncidv,varidv
    12741275             if (ncidpl.eq.-99) ncidpl=ncidv
    … …  
    12761277! Temperature
    12771278         if (guide_T) then
    1278              ncidt=NCOPN('T.nc',NCNOWRIT,rcode)
    1279              varidt=NCVID(ncidt,'AIR',rcode)
     1279             rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
     1280             rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
    12801281             print*,'ncidT,varidT',ncidt,varidt
    12811282             if (ncidpl.eq.-99) ncidpl=ncidt
    … …  
    12831284! Humidite
    12841285         if (guide_Q) then
    1285              ncidQ=NCOPN('hur.nc',NCNOWRIT,rcode)
    1286              varidQ=NCVID(ncidQ,'RH',rcode)
     1286             rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
     1287             rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
    12871288             print*,'ncidQ,varidQ',ncidQ,varidQ
    12881289             if (ncidpl.eq.-99) ncidpl=ncidQ
    … …  
    12901291! Pression de surface
    12911292         if ((guide_P).OR.(guide_modele)) then
    1292              ncidps=NCOPN('ps.nc',NCNOWRIT,rcode)
    1293              varidps=NCVID(ncidps,'SP',rcode)
     1293             rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
     1294             rcode = nf90_inq_varid(ncidps, 'SP', varidps)
    12941295             print*,'ncidps,varidps',ncidps,varidps
    12951296         endif
    12961297! Coordonnee verticale
    12971298         if (.not.guide_modele) then
    1298               rcode=NF_INQ_VARID(ncidpl,'LEVEL',varidpl)
    1299               IF (rcode.NE.0) rcode=NF_INQ_VARID(ncidpl,'PRESSURE',varidpl)
     1299              rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
     1300              IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
    13001301              print*,'ncidpl,varidpl',ncidpl,varidpl
    13011302         endif

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/ozonecm_m.F90

    r1182 r1188  
    1 !
    21! $Header$
    3 !
    4       SUBROUTINE ozonecm(rjour, rlat, paprs, o3)
    5       USE dimphy
    6       IMPLICIT none
    7 C
    8 C The ozone climatology is based on an analytic formula which fits the
    9 C Krueger and Mintzner (1976) profile, as well as the variations with
    10 C altitude and latitude of the maximum ozone concentrations and the total
    11 C column ozone concentration of Keating and Young (1986). The analytic
    12 C formula have been established by J-F Royer (CRNM, Meteo France), who
    13 C also provided us the code.
    14 C
    15 C A. J. Krueger and R. A. Minzner, A Mid-Latitude Ozone Model for the
    16 C 1976 U.S. Standard Atmosphere, J. Geophys. Res., 81, 4477, (1976).
    17 C
    18 C Keating, G. M. and D. F. Young, 1985: Interim reference models for the
    19 C middle atmosphere, Handbook for MAP, vol. 16, 205-229.
    20 C
     2module ozonecm_m
    213
    22 cym#include "dimensions.h"
    23 cym#include "dimphy.h"
    24 #include "clesphys.h"
    25 #include "YOMCST.h"
    26       REAL rlat(klon), paprs(klon,klev+1)
    27       REAL o3(klon,klev)   ! ozone concentration in kg/kg
    28       REAL tozon, rjour, pi, pl
    29       INTEGER i, k
    30 C----------------------------------------------------------
    31       REAL field(klon,klev+1)
    32       REAL ps
    33       PARAMETER (ps=101325.0)
    34       REAL an, unit, zo3q3
    35       SAVE an, unit, zo3q3
    36 c$OMP THREADPRIVATE(an, unit, zo3q3)
    37       REAL mu,gms, zslat, zsint, zcost, z, ppm, qpm, a
    38       REAL asec, bsec, aprim, zo3a3
    39 C----------------------------------------------------------
    40 c         data an /365.25/   (meteo)
    41       DATA an /360.00/
    42       DATA unit /2.1415e-05/
    43       DATA zo3q3 /4.0e-08/
     4  IMPLICIT NONE
    445
    45       pi = 4.0 * ATAN(1.0)
    46       DO k = 1, klev
    47       DO i = 1, klon
    48       zslat = SIN(pi/180.*rlat(i))
    49       zsint = SIN(2.*pi*(rjour+15.)/an)
    50       zcost = COS(2.*pi*(rjour+15.)/an)
    51       z = 0.0531+zsint*(-0.001595+0.009443*zslat) +
    52      .  zcost*(-0.001344-0.00346*zslat) +
    53      .  zslat**2*(.056222+zslat**2*(-.037609
    54      . +.012248*zsint+.00521*zcost+.008890*zslat))
    55       zo3a3 = zo3q3/ps/2.
    56       z = z-zo3q3*ps
    57       gms = z
    58       mu = ABS(sin(pi/180.*rlat(i)))
    59       ppm = 800.-(500.*zslat+150.*zcost)*zslat
    60       qpm = 1.74e-5-(7.5e-6*zslat+1.7e-6*zcost)*zslat
    61       bsec = 2650.+5000.*zslat**2
    62       a = 4.0*(bsec)**(3./2.)*(ppm)**(3./2.)*(1.0+(bsec/ps)**(3./2.))
    63       a = a/(bsec**(3./2.)+ppm**(3./2.))**2
    64       aprim = (2.666666*qpm*ppm-a*gms)/(1.0-a)
    65       aprim = amax1(0.,aprim)
    66       asec = (gms-aprim)*(1.0+(bsec/ps)**(3./2.))
    67       asec = AMAX1(0.0,asec)
    68       aprim = gms-asec/(1.+(bsec/ps)**(3./2.))
    69       pl = paprs(i,k)
    70       tozon = aprim/(1.+3.*(ppm/pl)**2)+asec/(1.+(bsec/pl)**(3./2.))
    71      .  + zo3a3*pl*pl
    72       tozon = tozon / 9.81  ! en kg/m**2
    73       tozon = tozon / unit  ! en kg/m**2  > u dobson (10e-5 m)
    74       tozon = tozon / 1000. ! en cm
    75       field(i,k) = tozon
    76       ENDDO
    77       ENDDO
    78 c
    79       DO i = 1, klon
    80          field(i,klev+1) = 0.0
    81       ENDDO
    82       DO k = 1, klev
    83         DO i = 1, klon
    84           o3(i,k) = field(i,k) - field(i,k+1)
    85         ENDDO
    86       ENDDO
    87 c
    88       RETURN
    89       END
     6contains
     7
     8  SUBROUTINE ozonecm(rjour,rlat,paprs,o3)
     9
     10    ! The ozone climatology is based on an analytic formula which fits the
     11    ! Krueger and Mintzner (1976) profile, as well as the variations with
     12    ! altitude and latitude of the maximum ozone concentrations and the total
     13    ! column ozone concentration of Keating and Young (1986). The analytic
     14    ! formula have been established by J.-F. Royer (CRNM, Meteo France), who
     15    ! also provided us the code.
     16
     17    ! A. J. Krueger and R. A. Minzner, A Mid-Latitude Ozone Model for the
     18    ! 1976 U.S. Standard Atmosphere, J. Geophys. Res., 81, 4477, (1976).
     19
     20    ! Keating, G. M. and D. F. Young, 1985: Interim reference models for the
     21    ! middle atmosphere, Handbook for MAP, vol. 16, 205-229.
     22
     23    USE dimphy, only: klon, klev
     24
     25    REAL, INTENT (IN) :: rjour
     26    REAL, INTENT (IN) :: rlat(klon), paprs(klon,klev+1)
     27    REAL o3(klon,klev)
     28    ! "o3(j, k)" is the column-density of ozone in cell "(j, k)", that is
     29    ! between interface "k" and interface "k + 1", in kDU.
     30
     31    ! Variables local to the procedure:
     32
     33    REAL tozon ! equivalent pressure of ozone above interface "k", in Pa
     34    real pi, pl
     35    INTEGER i, k
     36
     37    REAL field(klon,klev+1)
     38    ! "field(:, k)" is the column-density of ozone between interface
     39    ! "k" and the top of the atmosphere (interface "llm + 1"), in kDU.
     40
     41    real, PARAMETER:: ps=101325.
     42    REAL, parameter:: an = 360., zo3q3 = 4E-8
     43    REAL, parameter:: dobson_unit = 2.1415E-5 ! in kg m-2
     44    REAL gms, zslat, zsint, zcost, z, ppm, qpm, a
     45    REAL asec, bsec, aprim, zo3a3
     46
     47    !----------------------------------------------------------
     48
     49    pi = 4. * atan(1.)
     50    DO k = 1, klev
     51       DO i = 1, klon
     52          zslat = sin(pi / 180. * rlat(i))
     53          zsint = sin(2.*pi*(rjour+15.)/an)
     54          zcost = cos(2.*pi*(rjour+15.)/an)
     55          z = 0.0531 + zsint * (-0.001595+0.009443*zslat) &
     56               + zcost * (-0.001344-0.00346*zslat) &
     57               + zslat**2 * (.056222 + zslat**2 &
     58               * (-.037609+.012248*zsint+.00521*zcost+.008890*zslat))
     59          zo3a3 = zo3q3/ps/2.
     60          z = z - zo3q3*ps
     61          gms = z
     62          ppm = 800. - (500.*zslat+150.*zcost)*zslat
     63          qpm = 1.74E-5 - (7.5E-6*zslat+1.7E-6*zcost)*zslat
     64          bsec = 2650. + 5000.*zslat**2
     65          a = 4.0*(bsec)**(3./2.)*(ppm)**(3./2.)*(1.0+(bsec/ps)**(3./2.))
     66          a = a/(bsec**(3./2.)+ppm**(3./2.))**2
     67          aprim = (2.666666*qpm*ppm-a*gms)/(1.0-a)
     68          aprim = amax1(0., aprim)
     69          asec = (gms-aprim)*(1.0+(bsec/ps)**(3./2.))
     70          asec = amax1(0.0, asec)
     71          aprim = gms - asec/(1.+(bsec/ps)**(3./2.))
     72          pl = paprs(i, k)
     73          tozon = aprim / (1. + 3. * (ppm / pl)**2) &
     74               + asec / (1. + (bsec / pl)**(3./2.)) + zo3a3 * pl * pl
     75          ! Convert from Pa to kDU:
     76          field(i, k) = tozon / 9.81 / dobson_unit / 1e3
     77       END DO
     78    END DO
     79
     80    field(:,klev+1) = 0.
     81    forall (k = 1: klev) o3(:,k) = field(:,k) - field(:,k+1)
     82
     83  END SUBROUTINE ozonecm
     84
     85end module ozonecm_m

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F

    r1183 r1188  
    3434      use mod_phys_lmdz_mpi_data, only: is_mpi_root
    3535      USE aero_mod
     36      use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
    3637
    3738      IMPLICIT none
    … …  
    735736      EXTERNAL o3cm      ! initialiser l'ozone
    736737      EXTERNAL orbite    ! calculer l'orbite terrestre
    737       EXTERNAL ozonecm   ! prescrire l'ozone
    738738      EXTERNAL phyetat0  ! lire l'etat initial de la physique
    739739      EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique

  • LMDZ4/branches/LMDZ4-dev/output.def

    r1103 r1188  
    591591flag_ovap         =  2, 3, 4, 1, 1   
    592592name_ovap         =  ovap
    593 # Geopotential height
     593# ?
     594flag_ovapinit     =  2, 3, 4, 1, 1
     595name_ovapinit     =  ovapinit
     596# ?
    594597flag_wvapp        =  2, 10, 10, 10, 10   
    595598name_wvapp        =  wvapp
    596 # Zonal wind
     599# Geopotential height
    597600flag_geop         =  2, 3, 10, 1, 1   
    598601name_geop         =  geop
Note: See TracChangeset for help on using the changeset viewer.