Ignore:
Timestamp:
07/26/12 04:05:44 (12 years ago)
Author:
dubos
Message:

Minor changes :
caldyn_sw.f90, advect_tracer.f90
icosa_mod.f90 : added parameters for NCAR test cases needing global scope
guided_mod.f90 : CALL to guided_ncar now takes tt=it*dt instead of it as input

Significant changes :
timeloop_gcm.f90 : re-activated CALL to advection scheme
disvert_ncar.f90,
etat0_ncar.f90
guided_ncar_mod.f90 : simplification, introduced several getin(...), update due to recent changes in advection test cases (deformational flow, Hadley cell)
run_adv.def : new keys, reorganized for legibility

Tests :
icosa_gcm.exe tested with ncar_adv_shape=const and ncar_adv_wind=solid,deform,hadley.
q1=1 maintained to machine accuracy. Surface pressure slightly oscillates as expected.

FIXME : Tests by Sarvesh with revision 24 show incorrect advection of cosine bell by solid-body rotation. Not fixed.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/trunk/src/guided_ncar_mod.f90

    r19 r25  
    33  PRIVATE 
    44   
    5   REAL(rstd),SAVE :: dt 
    6   INTEGER,SAVE    :: test 
     5  INTEGER,SAVE :: case_wind 
    76   
     7  REAL(rstd), PARAMETER :: alpha=0.0 ! tilt of solid-body rotation 
     8  REAL(rstd), PARAMETER :: tau = 12*daysec ! 12 days               ! see p. 16 
     9  REAL(rstd), PARAMETER :: w0_deform = 23000*pi/tau, b=0.2, ptop=25494.4  ! see p. 16  
     10  REAL(rstd), PARAMETER :: u0_hadley=40.,w0_hadley=0.15 ,ztop= 12000.  
     11  INTEGER, PARAMETER    :: K_hadley=5 
     12 
    813  PUBLIC init_guided, guided 
    914 
    1015CONTAINS 
    1116 
    12  
    13   SUBROUTINE init_guided(dt_in) 
    14   IMPLICIT NONE 
    15     REAL(rstd),INTENT(IN) :: dt_in 
    16      
    17     dt=dt_in 
    18     test=2 
    19      
     17  SUBROUTINE init_guided 
     18    IMPLICIT NONE 
     19    CHARACTER(LEN=255) :: wind 
     20    wind='deform' 
     21    CALL getin('ncar_adv_wind',wind) 
     22    SELECT CASE(TRIM(wind)) 
     23    CASE('solid') 
     24       case_wind=0 
     25    CASE('deform') 
     26       case_wind=1 
     27    CASE('hadley') 
     28       case_wind=2 
     29    CASE DEFAULT 
     30       PRINT*,'Bad selector for variable ncar_adv_wind : <', TRIM(wind),'> options are <solid>, <deform>, <hadley>' 
     31       END SELECT 
    2032  END SUBROUTINE init_guided 
    21  
    2233   
    23   SUBROUTINE guided(it, f_ps, f_theta_rhodz, f_u, f_q) 
     34  SUBROUTINE guided(tt, f_ps, f_theta_rhodz, f_u, f_q) 
    2435  USE icosa 
    2536    IMPLICIT NONE 
    26     INTEGER, INTENT(IN)   :: it 
     37    REAL(rstd), INTENT(IN):: tt 
    2738    TYPE(t_field),POINTER :: f_ps(:) 
    2839    TYPE(t_field),POINTER :: f_phis(:) 
    … …  
    3849      CALL swap_geometry(ind)  
    3950      ue = f_u(ind)   
    40       CALL wind_profile(it,ue) 
     51      CALL wind_profile(tt,ue) 
    4152    END DO     
    4253 
    … …  
    4455   
    4556 
    46   SUBROUTINE wind_profile(it,ue) 
    47   USE icosa 
    48   IMPLICIT NONE 
    49     INTEGER,INTENT(IN) :: it  
     57  SUBROUTINE wind_profile(tt,ue) 
     58    USE icosa 
     59    USE disvert_mod 
     60    IMPLICIT NONE 
     61    REAL(rstd),INTENT(IN)  :: tt ! current time  
    5062    REAL(rstd),INTENT(OUT) :: ue(iim*3*jjm,llm) 
    5163    REAL(rstd) :: lon, lat 
    5264    REAL(rstd) :: nx(3),n_norm,Velocity(3,llm) 
    53     REAL(rstd), PARAMETER :: lon0=3*pi/2,lat0=0.0 
    5465    REAL(rstd) :: rr1,rr2,bb,cc,aa,hmx 
    5566    REAL(rstd) :: v1(3),v2(3),ny(3) 
    5667    INTEGER :: i,j,n,l 
    57     REAL(rstd) :: zrl(llm) 
     68    REAL(rstd) :: pitbytau,kk, pr, zr, u0, u1, v0 
    5869 
    59     CALL compute_zrl(zrl) 
     70    pitbytau = pi*tt/tau 
     71    kk = 10*radius/tau 
     72    u0 = 2*pi*radius/tau ! for solid-body rotation 
    6073!--------------------------------------------------------- 
    6174    DO l = 1,llm  
    62       DO j=jj_begin-1,jj_end+1 
    63         DO i=ii_begin-1,ii_end+1 
    64           n=(j-1)*iim+i 
    65           CALL compute_velocity(xyz_e(n+u_right,:),l,velocity(:,l)) 
    66           CALL cross_product2(xyz_v(n+z_rdown,:)/radius,xyz_v(n+z_rup,:)/radius,nx)        
    67           ue(n+u_right,l)=1e-10 
    68           n_norm=sqrt(sum(nx(:)**2)) 
    69           IF (n_norm>1e-30) THEN 
    70             nx=-nx/n_norm*ne(n,right) 
    71             ue(n+u_right,l)=sum(nx(:)*velocity(:,l)) 
    72             IF (ABS(ue(n+u_right,l))<1e-100) PRINT *,"ue(n+u_right) ==0",i,j,velocity(:,1) 
    73           ENDIF 
     75       pr = presnivs(l) 
     76       zr = -kappa*cpp*ncar_T0/g*log(pr/ncar_p0)  ! reciprocal of (1) p. 13, isothermal atmosphere 
     77       u1 = w0_deform*radius/b/ptop*cos(2*pitbytau)*(exp((ptop-pr)/b/ptop)-exp((pr-ncar_p0)/b/ptop)) 
     78       v0 = -radius*w0_hadley*pi/(5.0*ztop)*(ncar_p0/pr)*cos(pi*zr/ztop)*cos(pitbytau) ! for Hadley cell 
     79 
     80       DO j=jj_begin-1,jj_end+1 
     81          DO i=ii_begin-1,ii_end+1 
     82             n=(j-1)*iim+i 
     83             CALL compute_velocity(xyz_e(n+u_right,:),l,velocity(:,l)) 
     84             CALL cross_product2(xyz_v(n+z_rdown,:)/radius,xyz_v(n+z_rup,:)/radius,nx)        
     85             ue(n+u_right,l)=1e-10 
     86             n_norm=sqrt(sum(nx(:)**2)) 
     87             IF (n_norm>1e-30) THEN 
     88                nx=-nx/n_norm*ne(n,right) 
     89                ue(n+u_right,l)=sum(nx(:)*velocity(:,l)) 
     90                IF (ABS(ue(n+u_right,l))<1e-100) PRINT *,"ue(n+u_right)==0",i,j,velocity(:,1) 
     91             ENDIF 
    7492              
    75           CALL compute_velocity(xyz_e(n+u_lup,:),l,velocity(:,l)) 
    76           CALL cross_product2(xyz_v(n+z_up,:)/radius,xyz_v(n+z_lup,:)/radius,nx) 
    77  
    78           ue(n+u_lup,l)=1e-10 
    79           n_norm=sqrt(sum(nx(:)**2)) 
    80           IF (n_norm>1e-30) THEN 
    81             nx=-nx/n_norm*ne(n,lup) 
    82             ue(n+u_lup,l)=sum(nx(:)*velocity(:,l)) 
    83           ENDIF 
    84         
    85           CALL compute_velocity(xyz_e(n+u_ldown,:),l,velocity(:,l)) 
    86           CALL cross_product2(xyz_v(n+z_ldown,:)/radius,xyz_v(n+z_down,:)/radius,nx) 
    87  
    88           ue(n+u_ldown,l)=1e-10 
    89           n_norm=sqrt(sum(nx(:)**2)) 
    90           IF (n_norm>1e-30) THEN 
    91             nx=-nx/n_norm*ne(n,ldown) 
    92             ue(n+u_ldown,l)=sum(nx(:)*velocity(:,l)) 
    93             IF (ABS(ue(n+u_ldown,l))<1e-100) PRINT *,"ue(n+u_ldown) ==0",i,j 
    94           ENDIF         
    95         ENDDO 
    96       ENDDO 
     93             CALL compute_velocity(xyz_e(n+u_lup,:),l,velocity(:,l)) 
     94             CALL cross_product2(xyz_v(n+z_up,:)/radius,xyz_v(n+z_lup,:)/radius,nx) 
     95              
     96             ue(n+u_lup,l)=1e-10 
     97             n_norm=sqrt(sum(nx(:)**2)) 
     98             IF (n_norm>1e-30) THEN 
     99                nx=-nx/n_norm*ne(n,lup) 
     100                ue(n+u_lup,l)=sum(nx(:)*velocity(:,l)) 
     101             ENDIF 
     102              
     103             CALL compute_velocity(xyz_e(n+u_ldown,:),l,velocity(:,l)) 
     104             CALL cross_product2(xyz_v(n+z_ldown,:)/radius,xyz_v(n+z_down,:)/radius,nx) 
     105              
     106             ue(n+u_ldown,l)=1e-10 
     107             n_norm=sqrt(sum(nx(:)**2)) 
     108             IF (n_norm>1e-30) THEN 
     109                nx=-nx/n_norm*ne(n,ldown) 
     110                ue(n+u_ldown,l)=sum(nx(:)*velocity(:,l)) 
     111                IF (ABS(ue(n+u_ldown,l))<1e-100) PRINT *,"ue(n+u_ldown)==0",i,j 
     112             ENDIF 
     113          ENDDO 
     114       ENDDO 
    97115    END DO 
    98   
     116     
    99117  CONTAINS 
    100118            
    … …  
    104122      INTEGER,INTENT(IN)::l 
    105123      REAL(rstd),INTENT(OUT) :: velocity(3) 
    106       REAL(rstd) :: u0 
    107       REAL(rstd)  :: e_lat(3), e_lon(3) 
     124      REAL(rstd) :: e_lat(3), e_lon(3) 
    108125      REAL(rstd) :: lon,lat 
    109       REAL(rstd),PARAMETER::alpha=0.0 
    110126      REAL(rstd) :: u,v 
    111 !--------------------------------- test 1 
    112       REAL(rstd),parameter:: nday=12,daysec=86400 
    113       REAL(rstd),PARAMETER :: tau = nday*daysec 
    114       REAL(rstd) :: tt, pitbytau,kk 
    115 !---------------------------------- hadley  
    116       INTEGER,PARAMETER::K_hadly=5 
    117       REAL(rstd),PARAMETER::u0_hadly=40.,w0_hadly=0.25 ,ztop= 12000.  
    118       REAL(rstd)::coff_hadly 
    119          
    120          
    121       tt = it*dt  
    122       pitbytau = pi*tt/tau  
    123       kk = 10*radius/tau  
    124        
     127               
    125128      CALL xyz2lonlat(x/radius,lon,lat) 
    126129      e_lat(1) = -cos(lon)*sin(lat) 
    … …  
    133136 
    134137      u = 0.0 ; v = 0.0 
    135       u0 = 2*pi*radius/tau 
    136138     
    137       SELECT CASE(test)  
    138         CASE(0)   
    139           u=u0*(cos(lat)*cos(alpha)+sin(lat)*sin(alpha)*cos(lon)) 
    140           v=-u0*sin(lon)*sin(alpha) 
    141         CASE(1)  
    142           lon = lon - 2*pitbytau  
    143           u= kk*sin(lon)*sin(lon)*sin(2*lat)*cos(pitbytau)+ u0*cos(lat)  
    144           v= kk*sin(2*lon)*cos(lat)*cos(pitbytau)  
    145         CASE(2) 
    146           coff_hadly = -radius*w0_hadly*pi/(5.0*ztop)   
    147           u= u0_hadly*cos(lat)  
    148           v= coff_hadly*cos(lat)*sin(5.*lat)*cos(pi*zrl(l)/ztop)*cos(pitbytau) 
    149         CASE DEFAULT  
    150           PRINT*,"not valid choice of wind"  
    151        END SELECT  
    152  
     139      SELECT CASE(case_wind)  
     140      CASE(0)  ! Solid-body rotation 
     141         u=u0*(cos(lat)*cos(alpha)+sin(lat)*sin(alpha)*cos(lon)) 
     142         v=-u0*sin(lon)*sin(alpha) 
     143      CASE(1)  ! 3D Deformational flow -  
     144         lon = lon-2*pitbytau 
     145         u = kk*sin(lon)*sin(lon)*sin(2*lat)*cos(pitbytau)+ u0*cos(lat) 
     146         u = u + u1*cos(lon)*cos(lat)**2 
     147         v = kk*sin(2*lon)*cos(lat)*cos(pitbytau)  
     148      CASE(2) ! Hadley-like flow 
     149         u = u0_hadley*cos(lat) 
     150         v = v0*cos(lat)*sin(5.*lat)    ! Eq. 37 p. 19 
     151      CASE DEFAULT  
     152         PRINT*,"not valid choice of wind"  
     153      END SELECT 
     154       
    153155      Velocity(:)=(u*e_lon(:)+v*e_lat(:)+1e-50) 
    154  
     156       
    155157    END SUBROUTINE compute_velocity 
    156158     
    157     SUBROUTINE compute_zrl(zrl) 
    158     USE disvert_mod 
    159     IMPLICIT NONE 
    160       REAL(rstd),INTENT(OUT) :: zrl(llm) 
    161       INTEGER :: l 
    162       REAL(rstd) :: zr(llm+1) 
    163       REAl(rstd) :: pr 
    164       REAL(rstd), PARAMETER :: T0=300 
    165       REAL(rstd), PARAMETER :: R_d=287.0    
    166       REAL(rstd), PARAMETER :: ps0=100000.0 
    167  
    168        
    169       DO    l    = 1, llm+1 
    170         pr = ap(l) + bp(l)*ps0 
    171         zr(l) = -R_d*T0/g*log(pr/ps0)  
    172       ENDDO     
    173  
    174       DO l = 1, llm  
    175         zrl(l) = 0.5*(zr(l) + zr(l+1))  
    176       END DO 
    177      
    178     END SUBROUTINE compute_zrl 
    179    
    180159  END SUBROUTINE  wind_profile 
    181160 
Note: See TracChangeset for help on using the changeset viewer.