Ignore:
Timestamp:
Dec 22, 2009, 12:07:26 PM (15 years ago)
Author:
lguez
Message:

1) Replaced calls to "float" by calls to "real" and "dble" in
"phys_cosp". "call construct_cosp_gridbox(float(itap),..." was a bug
since the corresponding dummy argument has the type "double
precision". (And "float" is not in the Fortran standard.)

2) Modifications for the program "create_etat0_limit" only, "gcm" is
not touched:

2.1) Removed generic interface "startget" in module "startvar". Replaced
calls to "startget" by calls to "startget_phys2d", "startget_dyn" and
"startget_phys1d".

2.2) Simplified the interface of "startget_dyn" and made it more secure by
removing arguments for sizes of arrays and using assumed-shape array
arguments.

2.3) Corrected bug in call to "startget_dyn" for northward velocity. See
ticket #25 in Trac.

2.4) Collected procedures "inter_barxy", "inter_barx", "inter_bary",
"ord_coord" and "ord_coordm" in a module. "inter_barxy" is the only
public procedure in the module. Rewrote those five procedures: removed
arguments for sizes of arrays; used assumed-shape array arguments;
used array expressions; translated some spaghetti code in "inter_bary"
into structured code; corrected bug in "inter_bary", described by
ticket #26 in Trac.

2.5) Corrected a bug in "inter_bary": "y0" was initialized at 0
instead of -90. This bug made values at the south pole wrong. (This bug
has been acknowledeged by P. Le Van.)

2.6) Made the variable "champ" in procedure "limit_netcdf" an array of
rank 2. Thus it can be an argument in the call to the newly-secured
"inter_barxy".

2.7) The files "start.nc", "startphy.nc" and "limit.nc" are impacted
by this revision. There is a significant change in the variable "vcov"
of "start.nc". See
​http://web.lmd.jussieu.fr/~lglmd/Rev_1293/index.html. Note also the
changed value of "vcov" at the south pole. There is very little change
in other variables, in all three files.

Location:
LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3d
Files:
5 deleted
3 edited
1 copied

Legend:

Unmodified
Added
Removed

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3d/etat0_netcdf.F

    r1279 r1293  
    263263        varname = 'masque'
    264264        masque(:,:) = 0.0
    265         CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0,
    266      , jjm ,rlonu,rlatv , interbar )
     265        CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, masque,
     266     $       0.0, jjm ,rlonu,rlatv , interbar )
    267267        WRITE(*,*) 'MASQUE construit : Masque'
    268268        WRITE(*,'(97I1)') nINT(masque(:,:))
    … …  
    324324      ! This line needs to be replaced by a call to restget to get the values in the restart file
    325325      orog(:,:) = 0.0
    326        CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0 ,
    327      , jjm ,rlonu,rlatv , interbar, masque )
     326       CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, orog,
     327     $     0.0 , jjm ,rlonu,rlatv , interbar, masque )
    328328      !
    329329      WRITE(*,*) 'OUT OF GET VARIABLE : Relief'
    … …  
    333333      ! This line needs to be replaced by a call to restget to get the values in the restart file
    334334      rugo(:,:) = 0.0
    335        CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0 ,
    336      , jjm, rlonu,rlatv , interbar )
     335       CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, rugo,
     336     $     0.0 , jjm, rlonu,rlatv , interbar )
    337337      !
    338338      WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite'
    … …  
    346346      varname = 'psol'
    347347      psol(:,:) = 0.0
    348       CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0 ,
    349      , jjm ,rlonu,rlatv , interbar )
     348      CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, psol,
     349     $     0.0 , jjm ,rlonu,rlatv , interbar )
    350350      !
    351351      !  Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM
    … …  
    375375      varname = 'surfgeo'
    376376      phis(:,:) = 0.0
    377       CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0 ,
    378      , jjm ,rlonu,rlatv, interbar )
     377      CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, phis,
     378     $     0.0 , jjm ,rlonu,rlatv, interbar )
    379379      !
    380380      varname = 'u'
    381381      uvent(:,:,:) = 0.0
    382       CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls,
    383      . workvar, uvent, 0.0, jjm ,rlonv, rlatv, interbar )
     382      CALL startget_dyn(varname, rlonu, rlatu, pls, workvar, uvent, 0.,
     383     $     rlonv, rlatv, interbar )
    384384      ! 
    385385      varname = 'v'
    386386      vvent(:,:,:) = 0.0
    387       CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
    388      . workvar, vvent, 0.0, jjp1, rlonu, rlatu, interbar )
     387      CALL startget_dyn(varname, rlonv, rlatv, pls(:, :jjm, :),
     388     . workvar(:, :jjm, :), vvent, 0., rlonu, rlatu(:jjm), interbar )
    389389      !
    390390      varname = 't'
    391391      t3d(:,:,:) = 0.0
    392       CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
    393      . workvar, t3d, 0.0 , jjm, rlonu, rlatv , interbar )
     392      CALL startget_dyn(varname, rlonv, rlatu, pls, workvar, t3d, 0.,
     393     $    rlonu, rlatv , interbar )
    394394      !
    395395      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
    … …  
    397397      varname = 'tpot'
    398398      tpot(:,:,:) = 0.0
    399       CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
    400      . pk, tpot, 0.0 , jjm, rlonu, rlatv , interbar )
     399      CALL startget_dyn(varname, rlonv, rlatu, pls, pk, tpot, 0., rlonu,
     400     $     rlatv, interbar)
    401401      !
    402402      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
    … …  
    420420      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
    421421     .                           maxval(qsat(:,:,:))
    422       CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
    423      . qsat, qd, 0.0, jjm, rlonu, rlatv , interbar )
     422      CALL startget_dyn(varname, rlonv, rlatu, pls, qsat, qd, 0., rlonu,
     423     $    rlatv , interbar )
    424424      q3d(:,:,:,1) = qd(:,:,:)
    425425      !
    … …  
    431431      ! This line needs to be replaced by a call to restget to get the values in the restart file
    432432      tsol(:) = 0.0
    433       CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol, 0.0,
    434      .    jjm, rlonu, rlatv , interbar )
     433      CALL startget_phys1d(varname, iip1, jjp1, rlonv, rlatu, klon,
     434     $     tsol, 0.0, jjm, rlonu, rlatv , interbar )
    435435      !
    436436      WRITE(*,*) 'TSOL construit :'
    … …  
    439439      varname = 'qsol'
    440440      qsol(:) = 0.0
    441       CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol, 0.0,
    442      .  jjm, rlonu, rlatv , interbar )
     441      CALL startget_phys1d(varname, iip1, jjp1, rlonv, rlatu, klon,
     442     $     qsol, 0.0, jjm, rlonu, rlatv , interbar )
    443443      !
    444444      varname = 'snow'
    445445      sn(:) = 0.0
    446       CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn, 0.0,
    447      .    jjm, rlonu, rlatv , interbar )
     446      CALL startget_phys1d(varname, iip1, jjp1, rlonv, rlatu, klon, sn,
     447     $     0.0, jjm, rlonu, rlatv , interbar )
    448448      !
    449449      varname = 'rads'
    450450      radsol(:) = 0.0
    451       CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,
    452      .    jjm, rlonu, rlatv , interbar )
     451      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,
     452     $     0.0, jjm, rlonu, rlatv , interbar )
    453453      !
    454454      varname = 'rugmer'
    455455      rugmer(:) = 0.0
    456       CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,
    457      .    jjm, rlonu, rlatv , interbar )
     456      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,
     457     $     0.0, jjm, rlonu, rlatv , interbar )
    458458      !
    459459!      varname = 'agesno'
    460460!      agesno(:) = 0.0
    461 !      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0,
     461!      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0,
    462462!     .     jjm, rlonu, rlatv , interbar )
    463463
    464464      varname = 'zmea'
    465465      zmea(:) = 0.0
    466       CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,
     466      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,
    467467     .     jjm, rlonu, rlatv , interbar )
    468468
    469469      varname = 'zstd'
    470470      zstd(:) = 0.0
    471       CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,
     471      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,
    472472     .     jjm, rlonu, rlatv , interbar )
    473473      varname = 'zsig'
    474474      zsig(:) = 0.0
    475       CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,
     475      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,
    476476     .     jjm, rlonu, rlatv , interbar )
    477477      varname = 'zgam'
    478478      zgam(:) = 0.0
    479       CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,
     479      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,
    480480     .     jjm, rlonu, rlatv , interbar )
    481481      varname = 'zthe'
    482482      zthe(:) = 0.0
    483       CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,
     483      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,
    484484     .     jjm, rlonu, rlatv , interbar )
    485485      varname = 'zpic'
    486486      zpic(:) = 0.0
    487       CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,
     487      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,
    488488     .     jjm, rlonu, rlatv , interbar )
    489489      varname = 'zval'
    490490      zval(:) = 0.0
    491       CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,
     491      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,
    492492     .     jjm, rlonu, rlatv , interbar )
    493493c

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3d/inter_barxy_m.F90

    r1281 r1293  
    1 !
    2 ! $Header$
    3 !
    4        SUBROUTINE inter_barxy ( interfd,jnterfd,dlonid,dlatid ,
    5      ,        champ,imod,jmod,rlonimod,rlatimod, jsort,champint )
    6 
    7 c    Auteur :   P. Le Van
    8 c
    9        INTEGER interfd,jnterfd,imod,jmod
    10        REAL champ(interfd,jnterfd +1 ),dlonid(interfd),dlatid(jnterfd),
    11      ,      champint(imod,jsort)
    12        REAL rlonimod(imod),rlatimod(jmod)
    13 
    14 #include "dimensions.h"
    15 #include "paramet.h"
    16 #include "comgeom2.h"
    17 
    18        REAL champx(imod),champy(jnterfd +1,imod),chpn(imod),chps(imod)
    19        REAL chhpn,chhps
    20        REAL fmody(jjp1)
    21 c
    22 
    23        DO j = 1, jnterfd + 1
    24         CALL inter_barx( interfd, dlonid, champ( 1,j ),
    25      ,                       imod, rlonimod , champx )
    26          DO i = 1,imod
    27            champy(j,i) = champx(i)
    28          ENDDO
    29        ENDDO
    30 
    31        DO i = 1, imod
    32         CALL inter_bary( jjm,jnterfd,dlatid,champy(1,i),
    33      ,                     jmod ,rlatimod,  fmody     )
    34           DO j = 1, jsort
    35            champint(i,j) = fmody(j)
    36           ENDDO
    37        ENDDO
    38 
    39        IF( jsort.EQ.jjp1)  THEN
    40 
    41 c   ....  Valeurs uniques  aux  poles ....
    42 c
    43          DO i =  1,imod
    44           chpn(i)  = aire( i,  1   ) * champint( i, 1   )
    45           chps(i)  = aire( i, jjp1 ) * champint( i,jjp1 )
    46          ENDDO
    47           chhpn  = SSUM(imod,chpn,1)/apoln
    48           chhps  = SSUM(imod,chps,1)/apols
    49 
    50          DO i = 1, imod
    51           champint( i,  1  ) = chhpn
    52           champint( i, jjp1) = chhps
    53          ENDDO
    54 c
    55        ENDIF
    56 
    57        RETURN
    58        END
    59 
     1module inter_barxy_m
     2
     3  ! Authors: Robert SADOURNY, Phu LE VAN, Lionel GUEZ
     4
     5  implicit none
     6
     7  private
     8  public inter_barxy
     9
     10contains
     11
     12  SUBROUTINE inter_barxy(dlonid, dlatid, champ, rlonimod, rlatimod, champint)
     13
     14    use assert_eq_m, only: assert_eq
     15    use assert_m, only: assert
     16
     17    include "dimensions.h"
     18    ! (for "iim", "jjm")
     19
     20    include "paramet.h"
     21    ! (for other included files)
     22
     23    include "comgeom2.h"
     24    ! (for "aire", "apoln", "apols")
     25
     26    REAL, intent(in):: dlonid(:)
     27    ! (longitude from input file, in rad, from -pi to pi)
     28
     29    REAL, intent(in):: dlatid(:), champ(:, :), rlonimod(:)
     30
     31    REAL, intent(in):: rlatimod(:)
     32    ! (latitude angle, in degrees or rad, in strictly decreasing order)
     33
     34    real, intent(out):: champint(:, :)
     35    ! Si taille de la seconde dim = jjm + 1, on veut interpoler sur les
     36    ! jjm+1 latitudes rlatu du modele (latitudes des scalaires et de U)
     37    ! Si taille de la seconde dim = jjm, on veut interpoler sur les
     38    ! jjm latitudes rlatv du modèle (latitudes de V)
     39
     40    ! Variables local to the procedure:
     41
     42    REAL champy(iim, size(champ, 2))
     43    integer j, i, jnterfd, jmods
     44
     45    REAL yjmod(size(champint, 2))
     46    ! (angle, in degrees, in strictly increasing order)
     47
     48    REAL   yjdat(size(dlatid) + 1) ! angle, in degrees, in increasing order
     49    LOGICAL decrois ! "dlatid" is in decreasing order
     50
     51    !-----------------------------------
     52
     53    jnterfd = assert_eq(size(champ, 2) - 1, size(dlatid), &
     54         "inter_barxy jnterfd")
     55    jmods = size(champint, 2)
     56    call assert(size(champ, 1) == size(dlonid), "inter_barxy size(champ, 1)")
     57    call assert((/size(rlonimod), size(champint, 1)/) == iim, &
     58         "inter_barxy iim")
     59    call assert(any(jmods == (/jjm, jjm + 1/)), 'inter_barxy jmods')
     60    call assert(size(rlatimod) == jjm, "inter_barxy size(rlatimod)")
     61
     62    ! Check decreasing order for "rlatimod":
     63    DO i = 2, jjm
     64       IF (rlatimod(i) >= rlatimod(i-1)) stop &
     65            '"inter_barxy": "rlatimod" should be strictly decreasing'
     66    ENDDO
     67
     68    yjmod(:jjm) = ord_coordm(rlatimod)
     69    IF (jmods == jjm + 1) THEN
     70       IF (90. - yjmod(jjm) < 0.01) stop &
     71            '"inter_barxy": with jmods = jjm + 1, yjmod(jjm) should be < 90.'
     72    ELSE
     73       ! jmods = jjm
     74       IF (ABS(yjmod(jjm) - 90.) > 0.01) stop &
     75            '"inter_barxy": with jmods = jjm, yjmod(jjm) should be 90.'
     76    ENDIF
     77
     78    if (jmods == jjm + 1) yjmod(jjm + 1) = 90.
     79
     80    DO j = 1, jnterfd + 1
     81       champy(:, j) = inter_barx(dlonid, champ(:, j), rlonimod)
     82    ENDDO
     83
     84    CALL ord_coord(dlatid, yjdat, decrois)
     85    IF (decrois) champy(:, :) = champy(:, jnterfd + 1:1:-1)
     86    DO i = 1, iim
     87       champint(i, :) = inter_bary(yjdat, champy(i, :), yjmod)
     88    ENDDO
     89    champint(:, :) = champint(:, jmods:1:-1)
     90
     91    IF (jmods == jjm + 1) THEN
     92       ! Valeurs uniques aux poles
     93       champint(:, 1) = SUM(aire(:iim,  1) * champint(:, 1)) / apoln
     94       champint(:, jjm + 1) = SUM(aire(:iim, jjm + 1) &
     95            * champint(:, jjm + 1)) / apols
     96    ENDIF
     97
     98  END SUBROUTINE inter_barxy
     99
     100  !******************************
     101
     102  function inter_barx(dlonid, fdat, rlonimod)
     103
     104    !        INTERPOLATION BARYCENTRIQUE BASEE SUR LES AIRES
     105    !            VERSION UNIDIMENSIONNELLE  ,   EN  LONGITUDE .
     106
     107    !     idat : indice du champ de donnees, de 1 a idatmax
     108    !     imod : indice du champ du modele,  de 1 a  imodmax
     109    !     fdat(idat) : champ de donnees (entrees)
     110    !     inter_barx(imod) : champ du modele (sorties)
     111    !     dlonid(idat): abscisses des interfaces des mailles donnees
     112    !     rlonimod(imod): abscisses des interfaces des mailles modele
     113    !      ( L'indice 1 correspond a l'interface mailLE 1 / maille 2)
     114    !      ( Les abscisses sont exprimées en degres)
     115
     116    use assert_eq_m, only: assert_eq
     117
     118    IMPLICIT NONE
     119
     120    REAL, intent(in):: dlonid(:)
     121    real, intent(in):: fdat(:)
     122    real, intent(in):: rlonimod(:)
     123
     124    real inter_barx(size(rlonimod))
     125
     126    !    ...  Variables locales ...
     127
     128    INTEGER idatmax, imodmax
     129    REAL xxid(size(dlonid)+1), xxd(size(dlonid)+1), fdd(size(dlonid)+1)
     130    REAL  fxd(size(dlonid)+1), xchan(size(dlonid)+1), fdchan(size(dlonid)+1)
     131    REAL  xxim(size(rlonimod))
     132
     133    REAL x0, xim0, dx, dxm
     134    REAL chmin, chmax, pi
     135
     136    INTEGER imod, idat, i, ichang, id0, id1, nid, idatmax1
     137
     138    !-----------------------------------------------------
     139
     140    idatmax = assert_eq(size(dlonid), size(fdat), "inter_barx idatmax")
     141    imodmax = size(rlonimod)
     142
     143    pi = 2. * ASIN(1.)
     144
     145    !   REDEFINITION DE L'ORIGINE DES ABSCISSES
     146    !    A L'INTERFACE OUEST DE LA PREMIERE MAILLE DU MODELE 
     147    DO imod = 1, imodmax
     148       xxim(imod) = rlonimod(imod)
     149    ENDDO
     150
     151    CALL minmax( imodmax, xxim, chmin, chmax)
     152    IF( chmax.LT.6.50 )   THEN
     153       DO imod = 1, imodmax
     154          xxim(imod) = xxim(imod) * 180./pi
     155       ENDDO
     156    ENDIF
     157
     158    xim0 = xxim(imodmax) - 360.
     159
     160    DO imod = 1, imodmax
     161       xxim(imod) = xxim(imod) - xim0
     162    ENDDO
     163
     164    idatmax1 = idatmax +1
     165
     166    DO idat = 1, idatmax
     167       xxd(idat) = dlonid(idat)
     168    ENDDO
     169
     170    CALL minmax( idatmax, xxd, chmin, chmax)
     171    IF( chmax.LT.6.50 )  THEN
     172       DO idat = 1, idatmax
     173          xxd(idat) = xxd(idat) * 180./pi
     174       ENDDO
     175    ENDIF
     176
     177    DO idat = 1, idatmax
     178       xxd(idat) = AMOD( xxd(idat) - xim0, 360. )
     179       fdd(idat) = fdat (idat)
     180    ENDDO
     181
     182    i = 2
     183    DO while (xxd(i) >= xxd(i-1) .and. i < idatmax)
     184       i = i + 1
     185    ENDDO
     186    IF (xxd(i) < xxd(i-1)) THEN
     187       ichang = i
     188       !  ***  reorganisation  des longitudes entre 0. et 360. degres ****
     189       nid = idatmax - ichang +1
     190       DO i = 1, nid
     191          xchan (i) = xxd(i+ichang -1 )
     192          fdchan(i) = fdd(i+ichang -1 )
     193       ENDDO
     194       DO i=1, ichang -1
     195          xchan (i+ nid) = xxd(i)
     196          fdchan(i+nid) = fdd(i)
     197       ENDDO
     198       DO i =1, idatmax
     199          xxd(i) = xchan(i)
     200          fdd(i) = fdchan(i)
     201       ENDDO
     202    end IF
     203
     204    !    translation des champs de donnees par rapport
     205    !    a la nouvelle origine, avec redondance de la
     206    !       maille a cheval sur les bords
     207
     208    id0 = 0
     209    id1 = 0
     210
     211    DO idat = 1, idatmax
     212       IF ( xxd( idatmax1- idat ).LT.360.) exit
     213       id1 = id1 + 1
     214    ENDDO
     215
     216    DO idat = 1, idatmax
     217       IF (xxd(idat).GT.0.) exit
     218       id0 = id0 + 1
     219    END DO
     220
     221    IF( id1 /= 0 ) then
     222       DO idat = 1, id1
     223          xxid(idat) = xxd(idatmax - id1 + idat) - 360.
     224          fxd (idat) = fdd(idatmax - id1 + idat)     
     225       END DO
     226       DO idat = 1, idatmax - id1
     227          xxid(idat + id1) = xxd(idat)
     228          fxd (idat + id1) = fdd(idat)
     229       END DO
     230    end IF
     231
     232    IF(id0 /= 0) then
     233       DO idat = 1, idatmax - id0
     234          xxid(idat) = xxd(idat + id0)
     235          fxd (idat) = fdd(idat + id0)
     236       END DO
     237
     238       DO idat = 1, id0
     239          xxid (idatmax - id0 + idat) =  xxd(idat) + 360.
     240          fxd  (idatmax - id0 + idat) =  fdd(idat)   
     241       END DO
     242    else
     243       DO idat = 1, idatmax
     244          xxid(idat)  = xxd(idat)
     245          fxd (idat)  = fdd(idat)
     246       ENDDO
     247    end IF
     248    xxid(idatmax1) = xxid(1) + 360.
     249    fxd (idatmax1) = fxd(1)
     250
     251    !   initialisation du champ du modele
     252
     253    inter_barx(:) = 0.
     254
     255    ! iteration
     256
     257    x0   = xim0
     258    dxm  = 0.
     259    imod = 1
     260    idat = 1
     261
     262    do while (imod <= imodmax)
     263       do while (xxim(imod).GT.xxid(idat))
     264          dx   = xxid(idat) - x0
     265          dxm  = dxm + dx
     266          inter_barx(imod) = inter_barx(imod) + dx * fxd(idat)
     267          x0   = xxid(idat)
     268          idat = idat + 1
     269       end do
     270       IF (xxim(imod).LT.xxid(idat)) THEN
     271          dx   = xxim(imod) - x0
     272          dxm  = dxm + dx
     273          inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm
     274          x0   = xxim(imod)
     275          dxm  = 0.
     276          imod = imod + 1
     277       ELSE
     278          dx   = xxim(imod) - x0
     279          dxm  = dxm + dx
     280          inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm
     281          x0   = xxim(imod)
     282          dxm  = 0.
     283          imod = imod + 1
     284          idat = idat + 1
     285       END IF
     286    end do
     287
     288  END function inter_barx
     289
     290  !******************************
     291
     292  function inter_bary(yjdat, fdat, yjmod)
     293
     294    ! Interpolation barycentrique basée sur les aires.
     295    ! Version unidimensionnelle, en latitude.
     296    ! L'indice 1 correspond à l'interface maille 1 -- maille 2.
     297
     298    use assert_m, only: assert
     299
     300    IMPLICIT NONE
     301
     302    REAL, intent(in):: yjdat(:)
     303    ! (angles, ordonnées des interfaces des mailles des données, in
     304    ! degrees, in increasing order)
     305
     306    REAL, intent(in):: fdat(:) ! champ de données
     307
     308    REAL, intent(in):: yjmod(:)
     309    ! (ordonnées des interfaces des mailles du modèle)
     310    ! (in degrees, in strictly increasing order)
     311
     312    REAL inter_bary(size(yjmod)) ! champ du modèle
     313
     314    ! Variables local to the procedure:
     315
     316    REAL y0, dy, dym
     317    INTEGER jdat ! indice du champ de données
     318    integer jmod ! indice du champ du modèle
     319
     320    !------------------------------------
     321
     322    call assert(size(yjdat) == size(fdat), "inter_bary")
     323
     324    ! Initialisation des variables
     325    inter_bary(:) = 0.
     326    y0    = -90.
     327    dym   = 0.
     328    jmod  = 1
     329    jdat  = 1
     330
     331    do while (jmod <= size(yjmod))
     332       do while (yjmod(jmod) > yjdat(jdat))
     333          dy         = yjdat(jdat) - y0
     334          dym        = dym + dy
     335          inter_bary(jmod) = inter_bary(jmod) + dy * fdat(jdat)
     336          y0         = yjdat(jdat)
     337          jdat       = jdat + 1
     338       end do
     339       IF (yjmod(jmod) < yjdat(jdat)) THEN
     340          dy         = yjmod(jmod) - y0
     341          dym        = dym + dy
     342          inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym
     343          y0         = yjmod(jmod)
     344          dym        = 0.
     345          jmod       = jmod + 1
     346       ELSE
     347          ! {yjmod(jmod) == yjdat(jdat)}
     348          dy         = yjmod(jmod) - y0
     349          dym        = dym + dy
     350          inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym
     351          y0         = yjmod(jmod)
     352          dym        = 0.
     353          jmod       = jmod + 1
     354          jdat       = jdat + 1
     355       END IF
     356    end do
     357    ! Le test de fin suppose que l'interface 0 est commune aux deux
     358    ! grilles "yjdat" et "yjmod".
     359
     360  END function inter_bary
     361
     362  !******************************
     363
     364  SUBROUTINE ord_coord(xi, xo, decrois)
     365
     366    ! This procedure receives an array of latitudes.
     367    ! It converts them to degrees if they are in radians.
     368    ! If the input latitudes are in decreasing order, the procedure
     369    ! reverses their order.
     370    ! Finally, the procedure adds 90° as the last value of the array.
     371
     372    use assert_eq_m, only: assert_eq
     373
     374    IMPLICIT NONE
     375
     376    include "comconst.h"
     377    ! (for "pi")
     378
     379    REAL, intent(in):: xi(:)
     380    ! (latitude, in degrees or radians, in increasing or decreasing order)
     381    ! ("xi" should contain latitudes from pole to pole.
     382    ! "xi" should contain the latitudes of the boundaries of grid
     383    ! cells, not the centers of grid cells.
     384    ! So the extreme values should not be 90° and -90°.)
     385
     386    REAL, intent(out):: xo(:) ! angles in degrees
     387    LOGICAL, intent(out):: decrois
     388
     389    ! Variables  local to the procedure:
     390    INTEGER nmax, i
     391
     392    !--------------------
     393
     394    nmax = assert_eq(size(xi), size(xo) - 1, "ord_coord")
     395
     396    ! Check monotonicity:
     397    decrois = xi(2) < xi(1)
     398    DO i = 3, nmax
     399       IF (decrois .neqv. xi(i) < xi(i-1)) stop &
     400            '"ord_coord":  latitudes are not monotonic'
     401    ENDDO
     402
     403    IF (abs(xi(1)) < pi) then
     404       ! "xi" contains latitudes in radians
     405       xo(:nmax) = xi(:) * 180. / pi
     406    else
     407       ! "xi" contains latitudes in degrees
     408       xo(:nmax) = xi(:)
     409    end IF
     410
     411    IF (ABS(abs(xo(1)) - 90) < 0.001 .or. ABS(abs(xo(nmax)) - 90) < 0.001) THEN
     412       print *, "ord_coord"
     413       PRINT *, '"xi" should contain the latitudes of the boundaries of ' &
     414            // 'grid cells, not the centers of grid cells.'
     415       STOP
     416    ENDIF
     417
     418    IF (decrois) xo(:nmax) = xo(nmax:1:- 1)
     419    xo(nmax + 1) = 90.
     420
     421  END SUBROUTINE ord_coord
     422
     423  !***********************************
     424
     425  function ord_coordm(xi)
     426
     427    ! This procedure converts to degrees, if necessary, and inverts the
     428    ! order.
     429
     430    IMPLICIT NONE
     431
     432    include "comconst.h"
     433    ! (for "pi")
     434
     435    REAL, intent(in):: xi(:) ! angle, in rad or degrees
     436    REAL ord_coordm(size(xi)) ! angle, in degrees
     437
     438    !-----------------------------
     439
     440    IF (xi(1) < 6.5) THEN
     441       ! "xi" is in rad
     442       ord_coordm(:) = xi(size(xi):1:-1) * 180. / pi
     443    else
     444       ! "xi" is in degrees
     445       ord_coordm(:) = xi(size(xi):1:-1)
     446    ENDIF
     447
     448  END function ord_coordm
     449
     450end module inter_barxy_m

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3d/limit_netcdf.F

    r1279 r1293  
    99      USE dimphy
    1010      use phys_state_var_mod , ONLY : pctsrf
     11      use inter_barxy_m, only: inter_barxy
     12
    1113      IMPLICIT none
    1214c
    … …  
    7072      REAL , ALLOCATABLE :: dlon(:), dlat(:)
    7173      REAL , ALLOCATABLE :: dlon_ini(:), dlat_ini(:)
    72       REAL , ALLOCATABLE :: champ_msk(:), champ(:)
     74      REAL , ALLOCATABLE :: champ_msk(:), champ(:, :)
    7375      REAL , ALLOCATABLE :: work(:,:)
    7476
    … …  
    335337      ENDIF
    336338c
    337       ALLOCATE( champ(imdep*jmdep) )
     339      ALLOCATE( champ(imdep, jmdep) )
    338340
    339341      DO  200 l = 1, lmdep
    … …  
    364366
    365367       IF ( interbar )   THEN
    366          DO j = 1, imdep * jmdep
    367            champ(j) = LOG(champ(j))
    368          ENDDO
    369 
    370368        IF( l.EQ.1 )  THEN
    371369         WRITE(6,*) '-------------------------------------------------',
    … …  
    376374     ,'------------------------'
    377375        ENDIF
    378         CALL inter_barxy ( imdep,jmdep -1,dlon,dlat,champ ,
    379      ,                  iim,jjm,rlonu,rlatv, jjp1,champint )
     376        CALL inter_barxy(dlon, dlat(:jmdep -1), log(champ), rlonu(:iim),
     377     $       rlatv, champint)
    380378         DO j=1,jjp1
    381379          DO i=1,iim
    … …  
    556554      ENDIF
    557555c
    558       ALLOCATE ( champ(imdep*jmdep) )
     556      ALLOCATE ( champ(imdep, jmdep) )
    559557
    560558      DO l = 1, lmdep
    … …  
    604602cIM   ENDDO
    605603cIM22/02/2002
    606          CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
    607      ,     champ, iim, jjm, rlonu, rlatv, jjp1, champint )
     604       CALL inter_barxy (dlon, dlat(:jmdep -1), champ, rlonu(:iim),
     605     $      rlatv, champint)
    608606      ELSE
    609607         CALL sea_ice(imdep, jmdep, dlon, dlat, champ,
    … …  
    881879      ENDIF
    882880
    883        ALLOCATE( champ(imdep*jmdep) )
     881       ALLOCATE( champ(imdep, jmdep) )
    884882       IF( extrap )   THEN
    885883         ALLOCATE ( work(imdep,jmdep) )
    … …  
    924922     ,'------------------------'
    925923        ENDIF
    926        CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
    927      , champ, iim, jjm, rlonu, rlatv, jjp1, champint )
     924        CALL inter_barxy (dlon, dlat(:jmdep -1), champ, rlonu(:iim),
     925     $       rlatv, champint)
    928926      ELSE
    929927       CALL grille_m(imdep, jmdep, dlon, dlat, champ,
    … …  
    10821080      ENDIF
    10831081c
    1084       ALLOCATE ( champ(imdep*jmdep) )
     1082      ALLOCATE ( champ(imdep, jmdep) )
    10851083
    10861084      DO l = 1, lmdep
    … …  
    11201118        ENDIF
    11211119
    1122        CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
    1123      , champ, iim, jjm, rlonu, rlatv, jjp1, champint )
     1120       CALL inter_barxy(dlon, dlat(:jmdep-1), champ, rlonu(:iim),
     1121     $       rlatv, champint)
    11241122      ELSE
    11251123       CALL grille_m(imdep, jmdep, dlon, dlat, champ,

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3d/startvar.F

    r1279 r1293  
    99    !      There are three ways to access data from the database of atmospheric data which
    1010    !       can be used to initialize the model. This depends on the type of field which needs
    11     !       to be extracted. In any case the call should come after a restget and should be of the type :
    12     !                CALL startget(...)
    13     !
     11    !       to be extracted.
    1412    !       We will details the possible arguments to startget here :
    1513    !
    1614    !        - A 2D variable on the dynamical grid :
    17     !           CALL startget(varname, iml, jml, lon_in, lat_in, champ, val_ex, jml2, lon_in2, lat_in2, interbar )             
     15    !           CALL startget_phys2d(varname, iml, jml, lon_in, lat_in, champ, val_ex, jml2, lon_in2, lat_in2, interbar )             
    1816    !
    1917    !        - A 1D variable on the physical grid :
    20     !            CALL startget(varname, iml, jml, lon_in, lat_in, nbindex, champ, val_exp, jml2, lon_in2, lat_in2, interbar )
     18    !            CALL startget_phys1d(varname, iml, jml, lon_in, lat_in, nbindex, champ, val_exp, jml2, lon_in2, lat_in2, interbar )
    2119    !
    2220    !
    2321    !         - A 3D variable on the dynamical grid :
    24     !            CALL startget(varname, iml, jml, lon_in, lat_in, lml, pls, workvar, champ, val_exp, jml2, lon_in2, lat_in2, interbar )
     22    !            CALL startget_dyn(varname, iml, jml, lon_in, lat_in, lml, pls, workvar, champ, val_exp, jml2, lon_in2, lat_in2, interbar )
    2523    !
    2624    !
    … …  
    4644    !
    4745      PRIVATE
    48       PUBLIC startget
    49     !
    50     !
    51       INTERFACE startget
    52         MODULE PROCEDURE startget_phys2d, startget_phys1d, startget_dyn
    53       END INTERFACE
     46      public startget_phys2d, startget_phys1d, startget_dyn
    5447    !
    5548      INTEGER, SAVE :: fid_phys, fid_dyn
    … …  
    553546      SUBROUTINE start_init_phys( iml, jml, lon_in, lat_in, jml2,
    554547     .                 lon_in2, lat_in2 , interbar )
     548
     549      use inter_barxy_m, only: inter_barxy
    555550    !
    556551      INTEGER, INTENT(in) :: iml, jml ,jml2
    … …  
    648643        WRITE(6,*) '-------------------------------------------------',
    649644     ,'--------------'
    650         CALL inter_barxy ( iml_phys,jml_phys -1,lon_rad,lat_rad ,
    651      ,   var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var   )
     645        CALL inter_barxy(lon_rad, lat_rad(:jml_phys -1), var_ana,
     646     $       lon_in2(:iml-1), lat_in2(:jml-1), tmp_var)
    652647      ELSE
    653648        CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
    … …  
    674669        WRITE(6,*) '-------------------------------------------------',
    675670     ,'--------------'
    676         CALL inter_barxy ( iml_phys,jml_phys -1,lon_rad,lat_rad ,
    677      ,   var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var  )
     671        CALL inter_barxy(lon_rad, lat_rad(:jml_phys -1), var_ana,
     672     $       lon_in2(:iml-1), lat_in2(:jml-1), tmp_var)
    678673      ELSE
    679674        CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
    … …  
    691686    !
    692687    !
    693       SUBROUTINE startget_dyn(varname, iml, jml, lon_in, lat_in,
    694      . lml, pls, workvar, champ, val_exp,jml2, lon_in2, lat_in2 ,
     688      SUBROUTINE startget_dyn(varname, lon_in, lat_in,
     689     . pls, workvar, champ, val_exp, lon_in2, lat_in2 ,
    695690     ,  interbar )
     691
     692      use assert_eq_m, only: assert_eq
    696693    !
    697694    !   ARGUMENTS
    698695    !
    699       CHARACTER*(*), INTENT(in) :: varname
    700       INTEGER, INTENT(in) :: iml, jml, lml, jml2
    701       REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
    702       REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
    703       REAL, INTENT(in) :: pls(iml, jml, lml)
    704       REAL, INTENT(in) :: workvar(iml, jml, lml)
    705       REAL, INTENT(inout) :: champ(iml, jml, lml)
     696      CHARACTER(len=*), INTENT(in) :: varname
     697      REAL, INTENT(in) :: lon_in(:) ! dim(iml)
     698      REAL, INTENT(in) :: lat_in(:) ! dim(jml)
     699      REAL, INTENT(in) :: lon_in2(:) ! dim(iml)
     700      REAL, INTENT(in) :: lat_in2(:) ! dim(jml2)
     701      REAL, INTENT(in) :: pls(:, :, :) ! dim(iml, jml, lml)
     702      REAL, INTENT(in) :: workvar(:, :, :) ! dim(iml, jml, lml)
     703      REAL, INTENT(inout) :: champ(:, :, :) ! dim(iml, jml, lml)
    706704      REAL, INTENT(in) :: val_exp
    707705      LOGICAL interbar
    … …  
    709707    !    LOCAL
    710708    !
    711       INTEGER :: il, ij, ii
     709      INTEGER :: il, ij, ii, iml, jml, lml, jml2
    712710      REAL :: xppn, xpps
    713711    !
    … …  
    719717    !   This routine only works if the variable does not exist or is constant
    720718    !
     719C     -----------------------------
     720
     721      iml = assert_eq((/size(lon_in), size(pls, 1), size(workvar, 1),
     722     $     size(champ, 1), size(lon_in2)/), "startget_dyn iml")
     723      jml = assert_eq(size(lat_in), size(pls, 2), size(workvar, 2),
     724     $     size(champ, 2), "startget_dyn jml")
     725      lml = assert_eq(size(pls, 3), size(workvar, 3), size(champ, 3),
     726     $     "startget_dyn lml")
     727      jml2 = size(lat_in2)
     728
    721729      IF ( MINVAL(champ(:,:,:)).EQ.MAXVAL(champ(:,:,:)) .AND.
    722730     . MINVAL(champ(:,:,:)).EQ.val_exp ) THEN
    … …  
    834842     ,             lat_in2 , interbar )
    835843    !
     844      use inter_barxy_m, only: inter_barxy
     845
    836846      INTEGER, INTENT(in) :: iml, jml, jml2
    837847      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
    … …  
    929939        WRITE(6,*) '-------------------------------------------------',
    930940     ,'--------------'
    931         CALL inter_barxy ( iml_dyn,jml_dyn -1,lon_rad,lat_rad ,
    932      ,    var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var)
     941        CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1), var_ana,
     942     $       lon_in2(:iml-1), lat_in2(:jml-1), tmp_var)
    933943      ELSE
    934944        CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
    … …  
    954964        WRITE(6,*) '-------------------------------------------------',
    955965     ,'--------------'
    956         CALL inter_barxy ( iml_dyn,jml_dyn -1,lon_rad,lat_rad ,
    957      ,    var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var)
     966        CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1), var_ana,
     967     $       lon_in2(:iml-1), lat_in2(:jml-1), tmp_var)
    958968      ELSE
    959969        CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
    … …  
    10351045    !
    10361046    !
     1047      use inter_barxy_m, only: inter_barxy
     1048
    10371049    !    ARGUMENTS
    10381050    !
    … …  
    11201132     ,'--------------'
    11211133       ENDIF
    1122        CALL inter_barxy ( iml_dyn, jml_dyn -1,lon_rad, lat_rad,
    1123      , var_ana3d(:,:,il),iml-1, jml2, lon_in2, lat_in2,jml,var_tmp2d )
     1134       CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1),
     1135     $      var_ana3d(:,:,il), lon_in2(:iml-1), lat_in2, var_tmp2d)
    11241136      ELSE
    11251137       CALL grille_m(iml_dyn, jml_dyn, lon_rad, lat_rad,
Note: See TracChangeset for help on using the changeset viewer.