Ignore:
Timestamp:
Aug 22, 2013, 4:02:07 PM (11 years ago)
Author:
emillour
Message:

Common dynamics: (and collateral adaptations in Venus physics)
Improved cpdet routines in and additional sponge mode:

  • Additionnal sponge mode (trigered with "callsponge" flag), in line with the one used in the Generic and Martian GCM. This sponge is called whenever there is a dissipation step.
  • Improvement of the cpdet routines : created routines tpot2t_glo_p and t2tpot_glo_p which handle fields on the whole dynamics (scaler) grid, which are more efficient than calling tpot2t_p or t2tpot_p with slabs of data (generated use of intermediate copies of these chunks of data at every call)
  • Turned cpdet.F into a module cpdet_mod.F90 (and correspondingly adapted all routines in the Venus physics).

EM

Location:
trunk/LMDZ.COMMON/libf/dyn3dpar
Files:
1 added
5 edited
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F

    r892 r1017  
    3939      USE infotrac
    4040      USE control_mod
     41      USE cpdet_mod, only: tpot2t_p, t2tpot_p
    4142
    4243      IMPLICIT NONE

  • trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F

    r1012 r1017  
    1919      USE infotrac, ONLY : type_trac
    2020      use assert_m, only: assert
     21      use sponge_mod_p, only: callsponge,mode_sponge,nsponge,tetasponge
    2122      IMPLICIT NONE
    2223c-----------------------------------------------------------------------
    … …  
    388389       tau_top_bound=1.e-5
    389390       CALL getin('tau_top_bound',tau_top_bound)
     391
     392! the other possible sponge layer (sponge_mod)
     393       callsponge=.false. ! default value; don't use the sponge
     394       call getin("callsponge",callsponge)
     395       ! check that user is not trying to use both sponge models
     396       if ((iflag_top_bound.ge.1).and.callsponge) then
     397         write(lunout,*)'Bad choice of options:'
     398         write(lunout,*)' iflag_top_bound=',iflag_top_bound
     399         write(lunout,*)' and callsponge=.true.'
     400         write(lunout,*)'But both sponge models should not be',
     401     &                  ' used simultaneously!'
     402         stop
     403       endif
     404       
     405! nsponge: number of atmospheric layers over which the sponge extends
     406       nsponge=3 ! default value
     407       call getin("nsponge",nsponge)
     408
     409! mode_sponge: (quenching is towards ... over the upper nsponge layers)
     410!      0: (h=hmean,u=v=0)
     411!      1: (h=hmean,u=umean,v=0)
     412!      2: (h=hmean,u=umean,v=vmean)"
     413       mode_sponge=2 ! default value
     414       call getin("mode_sponge",mode_sponge)
     415
     416! tetasponge: characteristic time scale (seconds) at topmost layer
     417!            (time scale then doubles with decreasing layer index)."
     418       tetasponge=50000.0
     419       call getin("tetasponge",tetasponge)
    390420
    391421! FOR TITAN: tidal forces

  • trunk/LMDZ.COMMON/libf/dyn3dpar/cpdet_mod.F90

    r1016 r1017  
    1 ! ADAPTATION GCM POUR CP(T)
    2 c======================================================================
    3 c S. Lebonnois, 10/2010
    4 c
    5 c Cp doit être calculé par cpdet(t) pour être valable partout
    6 c
    7 c La fonction d'Exner reste pk = RCPD*(play/pref)**RKAPPA
    8 c (RCPD=cpp, RKAPPA=kappa)
    9 c
    10 c On passe de T a teta (temperature potentielle) par t2tpot(t,teta,pk)
    11 c On passe de teta a T par tpot2t(teta,t,pk)
    12 c
    13 c======================================================================
     1module cpdet_mod
     2
     3implicit none
     4
     5! ADAPTATION OF GCM TO CP(T)
     6!======================================================================
     7! S. Lebonnois, 10/2010
     8!
     9! Cp must be computed using cpdet(t) to be valid
     10!
     11! The Exner function is still pk = RCPD*(play/pref)**RKAPPA
     12! (RCPD=cpp, RKAPPA=kappa)
     13!
     14! One goes from T to teta (potential temperature) using t2tpot(t,teta,pk)
     15! One goes from teta to T using tpot2t(teta,t,pk)
     16!
     17!======================================================================
     18
     19contains
    1420
    1521      SUBROUTINE ini_cpdet
    … …  
    1723      USE control_mod, ONLY: planet_type
    1824      IMPLICIT none
    19 c======================================================================
    20 c Initialisation de nu_venus et t0_venus
    21 c======================================================================
     25!======================================================================
     26! Initialization of nu_venus and t0_venus
     27!======================================================================
    2228
    2329! for cpp, nu_venus and t0_venus:
    … …  
    3339
    3440      return
    35       end
    36 
    37 c======================================================================
    38 c======================================================================
     41      end subroutine ini_cpdet
     42
     43!======================================================================
     44!======================================================================
    3945
    4046      FUNCTION cpdet(t)
    … …  
    4652#include "comconst.h"
    4753
    48       real cpdet,t
     54      real,intent(in) :: t
     55      real cpdet
    4956
    5057      if (planet_type.eq."venus") then
    … …  
    5562
    5663      return
    57       end
    58      
    59 c======================================================================
    60 c======================================================================
     64      end function cpdet
     65     
     66!======================================================================
     67!======================================================================
    6168
    6269      SUBROUTINE t2tpot(npoints, yt, yteta, ypk)
    63 c======================================================================
    64 c Arguments:
    65 c
    66 c yt   --------input-R- Temperature
    67 c yteta-------output-R- Temperature potentielle
    68 c ypk  --------input-R- Fonction d'Exner: RCPD*(pplay/pref)**RKAPPA
    69 c
    70 c======================================================================
     70!======================================================================
     71! Arguments:
     72!
     73! yt   --------input-R- Temperature
     74! yteta-------output-R- Temperature potentielle
     75! ypk  --------input-R- Fonction d'Exner: RCPD*(pplay/pref)**RKAPPA
     76!
     77!======================================================================
    7178
    7279      USE control_mod, ONLY: planet_type
    … …  
    7683#include "comconst.h"
    7784
    78       integer npoints
    79       REAL    yt(npoints), yteta(npoints), ypk(npoints)
     85      integer,intent(in) :: npoints
     86      REAL,intent(in) :: yt(npoints), ypk(npoints)
     87      REAL,intent(out) :: yteta(npoints)
    8088     
    8189      if (planet_type.eq."venus") then
    … …  
    8896
    8997      return
    90       end
    91 
    92 c======================================================================
    93 c======================================================================
     98      end subroutine t2tpot
     99
     100!======================================================================
     101!======================================================================
    94102
    95103      SUBROUTINE t2tpot_p(nlon,nlev, yt, yteta, ypk)
    96 ! Parallel version of t2tpot
    97       USE parallel
     104! Parallel version of t2tpot, for an arbitrary number of columns
    98105      USE control_mod, only : planet_type
    99106      IMPLICIT none
    … …  
    110117     
    111118      if (planet_type.eq."venus") then
     119!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    112120        do l=1,nlev
    113121          yteta(:,l)=yt(:,l)**nu_venus                                  &
    … …  
    116124          yteta(:,l)=yteta(:,l)**(1./nu_venus)
    117125        enddo
    118       else
     126!$OMP END DO
     127      else
     128!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    119129        do l=1,nlev
    120130          yteta(:,l)=yt(:,l)*cpp/ypk(:,l)
    121131        enddo
     132!$OMP END DO
    122133      endif ! of if (planet_type.eq."venus")
    123134
    124       END
    125 
    126 c======================================================================
     135      end subroutine t2tpot_p
     136
     137!======================================================================
     138!======================================================================
     139
     140      SUBROUTINE t2tpot_glo_p(yt, yteta, ypk)
     141! Parallel version of t2tpot, over the full dynamics (scalar) grid
     142! (more efficient than multiple calls to t2tpot_p() with slices of data)
     143      USE parallel, only : jj_begin,jj_end
     144      USE control_mod, only : planet_type
     145      IMPLICIT none
     146! for iip1, jjp1 and llm
     147#include "dimensions.h"
     148#include "paramet.h"
     149! for cpp, nu_venus and t0_venus:
     150#include "comconst.h"
     151
     152      real,intent(in) :: yt(iip1,jjp1,llm)
     153      real,intent(out) :: yteta(iip1,jjp1,llm)
     154      real,intent(in) :: ypk(iip1,jjp1,llm)
     155! local variable:
     156      integer :: j,l
     157      integer :: jjb,jje
     158     
     159      jjb=jj_begin
     160      jje=jj_end
     161
     162      if (planet_type.eq."venus") then
     163!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     164        do l=1,llm
     165          yteta(:,jjb:jje,l)=yt(:,jjb:jje,l)**nu_venus                  &
     166     &                     -nu_venus*t0_venus**nu_venus*                &
     167     &                          log(ypk(:,jjb:jje,l)/cpp)
     168          yteta(:,jjb:jje,l)=yteta(:,jjb:jje,l)**(1./nu_venus)
     169        enddo
     170!$OMP END DO
     171      else
     172!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     173        do l=1,llm
     174          yteta(:,jjb:jje,l)=yt(:,jjb:jje,l)*cpp/ypk(:,jjb:jje,l)
     175        enddo
     176!$OMP END DO
     177      endif ! of if (planet_type.eq."venus")
     178
     179      end subroutine t2tpot_glo_p
     180
     181!======================================================================
    127182
    128183      SUBROUTINE tpot2t(npoints,yteta, yt, ypk)
    129 c======================================================================
    130 c Arguments:
    131 c
    132 c yteta--------input-R- Temperature potentielle
    133 c yt   -------output-R- Temperature
    134 c ypk  --------input-R- Fonction d'Exner: RCPD*(pplay/pref)**RKAPPA
    135 c
    136 c======================================================================
     184!======================================================================
     185! Arguments:
     186!
     187! yteta--------input-R- Temperature potentielle
     188! yt   -------output-R- Temperature
     189! ypk  --------input-R- Fonction d'Exner: RCPD*(pplay/pref)**RKAPPA
     190!
     191!======================================================================
    137192
    138193      USE control_mod, ONLY: planet_type
    … …  
    142197#include "comconst.h"
    143198
    144       integer npoints
    145       REAL    yt(npoints), yteta(npoints), ypk(npoints)
     199      integer,intent(in) :: npoints
     200      REAL,intent(in) :: yteta(npoints), ypk(npoints)
     201      REAL,intent(out) :: yt(npoints)
    146202     
    147203      if (planet_type.eq."venus") then
    … …  
    154210 
    155211      return
    156       end
    157 
    158 c======================================================================
    159 c======================================================================
     212      end subroutine tpot2t
     213
     214!======================================================================
     215!======================================================================
     216
    160217      SUBROUTINE tpot2t_p(nlon,nlev,yteta,yt,ypk)
    161 ! Parallel version of tpot2t
    162       USE parallel
     218! Parallel version of tpot2t, for an arbitrary number of columns
    163219      USE control_mod, only : planet_type
    164220      IMPLICIT none
    … …  
    175231
    176232      if (planet_type.eq."venus") then
     233!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    177234        do l=1,nlev
    178235          yt(:,l)=yteta(:,l)**nu_venus                                  &
    … …  
    181238          yt(:,l)=yt(:,l)**(1./nu_venus)
    182239        enddo
    183       else
     240!$OMP END DO
     241      else
     242!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    184243        do l=1,nlev
    185244          yt(:,l)=yteta(:,l)*ypk(:,l)/cpp
    186245        enddo
     246!$OMP END DO
    187247      endif ! of if (planet_type.eq."venus")
    188       END
    189 
    190 c======================================================================
    191 c======================================================================
    192 c
    193 c ATTENTION
    194 c
    195 c Si un jour on a besoin, il faudra coder les routines
    196 c    dt2dtpot / dtpto2dt
    197 c
    198 c======================================================================
    199 c======================================================================
     248      end subroutine tpot2t_p
     249
     250!======================================================================
     251!======================================================================
     252
     253      SUBROUTINE tpot2t_glo_p(yteta,yt,ypk)
     254! Parallel version of tpot2t, over the full dynamics (scalar) grid
     255! (more efficient than multiple calls to tpot2t_p() with slices of data)
     256      USE parallel, only : jj_begin,jj_end
     257      USE control_mod, only : planet_type
     258      IMPLICIT none
     259! for iip1, jjp1 and llm
     260#include "dimensions.h"
     261#include "paramet.h"
     262! for cpp, nu_venus and t0_venus:
     263#include "comconst.h"
     264
     265      real,intent(out) :: yt(iip1,jjp1,llm)
     266      real,intent(in) :: yteta(iip1,jjp1,llm)
     267      real,intent(in) :: ypk(iip1,jjp1,llm)
     268! local variable:
     269      integer :: j,l
     270      integer :: jjb,jje
     271     
     272      jjb=jj_begin
     273      jje=jj_end
     274
     275      if (planet_type.eq."venus") then
     276!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     277        do l=1,llm
     278          yt(:,jjb:jje,l)=yteta(:,jjb:jje,l)**nu_venus                  &
     279     &                  +nu_venus*t0_venus**nu_venus*                   &
     280     &                       log(ypk(:,jjb:jje,l)/cpp)
     281          yt(:,jjb:jje,l)=yt(:,jjb:jje,l)**(1./nu_venus)
     282        enddo
     283!$OMP END DO
     284      else
     285!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     286        do l=1,llm
     287          yt(:,jjb:jje,l)=yteta(:,jjb:jje,l)*ypk(:,jjb:jje,l)/cpp
     288        enddo
     289!$OMP END DO
     290      endif ! of if (planet_type.eq."venus")
     291      end subroutine tpot2t_glo_p
     292
     293!======================================================================
     294!======================================================================
     295!
     296! ATTENTION
     297!
     298! Si un jour on a besoin, il faudra coder les routines
     299!    dt2dtpot / dtpto2dt
     300!
     301!======================================================================
     302!======================================================================
     303end module cpdet_mod

  • trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F

    r979 r1017  
    1919      USE filtreg_mod
    2020      USE control_mod
     21      use cpdet_mod, only: ini_cpdet
    2122
    2223! Ehouarn: the following are needed with (parallel) physics:

  • trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F

    r1012 r1017  
    2121       USE getparam
    2222       USE control_mod
     23       use cpdet_mod, only: cpdet,tpot2t_glo_p,t2tpot_glo_p
     24       use sponge_mod_p, only: callsponge,mode_sponge,sponge_p
    2325
    2426      IMPLICIT NONE
    … …  
    186188      ! for CP(T)  -- Aymeric
    187189      real :: dtec
    188       real,external :: cpdet
    189190      real,save :: ztetaec(ip1jmp1,llm)  !!SAVE ???
    190191
    … …  
    211212c dummy: sinon cette routine n'est jamais compilee...
    212213      if(1.eq.0) then
     214#ifdef CPP_PHYS
    213215        CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
     216#endif
    214217      endif
    215218
    … …  
    550553c   --------------------------------
    551554! ADAPTATION GCM POUR CP(T)
     555      call tpot2t_glo_p(teta,temp,pk)
    552556      ijb=ij_begin
    553557      ije=ij_end
    554       call tpot2t_p(ije-ijb+1,llm,teta(ijb:ije,:),temp(ijb:ije,:),
    555      &                            pk(ijb:ije,:))
    556558!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    557559      do l=1,llm
    … …  
    11581160
    11591161      IF(apdiss) THEN
    1160 cc$OMP  PARALLEL DEFAULT(SHARED)
    1161 cc$OMP+          PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
     1162
    11621163c$OMP MASTER
    11631164        call suspend_timer(timer_caldyn)
    … …  
    11691170        call VTb(VThallo)
    11701171c$OMP END MASTER
     1172
     1173        ! sponge layer
     1174        if (callsponge) then
     1175          call Register_SwapFieldHallo(ps,ps,ip1jmp1,1,
     1176     &                          jj_Nb_dissip,1,1,Request_dissip)
     1177          call Register_Hallo(ps,ip1jm,1,1,1,1,1,Request_Dissip)
     1178          call SendRequest(Request_Dissip)
     1179c$OMP BARRIER
     1180          call WaitRequest(Request_Dissip)
     1181c$OMP BARRIER
     1182c$OMP MASTER
     1183          call VTe(VThallo)
     1184          call VTb(VThallo)
     1185c$OMP END MASTER
     1186c$OMP BARRIER
     1187          CALL sponge_p(ucov,vcov,teta,ps,dtdiss,mode_sponge)
     1188        endif
     1189
    11711190
    11721191c$OMP BARRIER
    … …  
    12051224c   dissipation
    12061225! ADAPTATION GCM POUR CP(T)
    1207             ijb=ij_begin
    1208             ije=ij_end
    1209         call tpot2t_p(ije-ijb+1,llm,teta(ijb:ije,:),temp(ijb:ije,:),
    1210      &                            pk(ijb:ije,:))
     1226        call tpot2t_glo_p(teta,temp,pk)
    12111227
    12121228!        CALL FTRACE_REGION_BEGIN("dissip")
    … …  
    12661282            enddo
    12671283c$OMP END DO
    1268         call t2tpot_p(ije-ijb+1,llm,temp(ijb:ije,:),ztetaec(ijb:ije,:),
    1269      &                            pk(ijb:ije,:))
     1284!        call t2tpot_p(ije-ijb+1,llm,temp(ijb:ije,:),ztetaec(ijb:ije,:),
     1285!     &                            pk(ijb:ije,:))
     1286         call t2tpot_glo_p(temp,ztetaec,pk)
    12701287c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)           
    12711288            do l=1,llm
    … …  
    15761593
    15771594! ADAPTATION GCM POUR CP(T)
     1595      call tpot2t_glo_p(teta,temp,pk)
    15781596      ijb=ij_begin
    15791597      ije=ij_end
    1580       call tpot2t_p(ije-ijb+1,llm,teta(ijb:ije,:),temp(ijb:ije,:),
    1581      &                            pk(ijb:ije,:))
    15821598!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    15831599      do l=1,llm
    … …  
    18011817
    18021818! ADAPTATION GCM POUR CP(T)
     1819                call tpot2t_glo_p(teta,temp,pk)
    18031820                ijb=ij_begin
    18041821                ije=ij_end
    1805            call tpot2t_p(ije-ijb+1,llm,teta(ijb:ije,:),temp(ijb:ije,:),
    1806      &                            pk(ijb:ije,:))
    18071822!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    18081823                do l=1,llm     

  • trunk/LMDZ.COMMON/libf/dyn3dpar/vlspltqs_p.F

    r109 r1017  
    2525      USE mod_hallo
    2626      USE VAMPIR
     27      use cpdet_mod, only: tpot2t_glo_p
    2728      IMPLICIT NONE
    2829
    … …  
    9798! ADAPTATION GCM POUR CP(T)
    9899! probablement a revoir...
    99       call tpot2t_p(ip1jmp1,llm,teta,tempe,pk)
    100 
     100!      call tpot2t_p(ip1jmp1,llm,teta,tempe,pk)
     101      call tpot2t_glo_p(teta,tempe,pk)
    101102
    102103      call SetTag(MyRequest1,100)
Note: See TracChangeset for help on using the changeset viewer.