Changeset 25


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.

Location:
codes/icosagcm/trunk
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/trunk/run_adv.def

    r18 r25  
     1#------------------------------- Mesh --------------------------------- 
     2 
    13# Number of subdivision on a main triangle (nbp) : integer (default=40) 
    24nbp=20 
    35 
    4 # Number of vertical layer (llm) : integer (default=19) 
     6# optim_it : mesh optimisation : number of iteration : integer (default=0) 
     7optim_it=0 
     8 
     9# Number of vertical layers (llm) : integer (default=19) 
    510llm=19 
     11 
     12# disvert : vertical discretisation : string (default='std') : std, ncar 
     13disvert=ncar 
     14 
     15#---------------------------------- Misc -------------------------------- 
    616 
    717# number of tracer (nqtot) : integer (default 1) 
    … …  
    1121dt = 180. 
    1222 
    13 # scheme type : string ( default='adam_bashforth') euler, leapfrog, leapfrog_matsuno, adam_bashforth) 
     23# number of timestep (default 100) 
     24itaumax = 38400 
     25 
     26# output field period : integer (default none) 
     27write_period=3600 
     28 
     29# etat0 : initial state : string (default=jablonowsky06) :  
     30# jablonowsky06, academic, ncar 
     31etat0=ncar 
     32 
     33# caldyn : computation type for gcm equation : string (default=gcm) : gcm, adv 
     34caldyn=adv 
     35 
     36#------------------------------ Dynamics -------------------------------- 
     37 
     38# scheme type : string ( default='adam_bashforth') 
     39# euler, leapfrog, leapfrog_matsuno, adam_bashforth) 
    1440scheme = euler 
    1541 
    1642# matsuno period : integer ( default=5) 
    1743matsuno_period = 5 
    18  
    19 # number of timestep (default 100) 
    20 itaumax = 38400 
    2144 
    2245# dissipation time graddiv : real (default=5000) 
    … …  
    2548 
    2649# dissipation time nxgradrot (default=5000) 
    27 tetagrot = 1800 
    28 tau_gradrot = 1800 
    29  
     50tetagrot=1800 
     51tau_gradrot=1800 
    3052tau_divgrad=1800 
    31  
    32 # output field period : integer (default none) 
    33 write_period=3600 
    34  
    35  
    36 #NCAR testcase : integer ( default=5) 
    37 ncar_test_case=1 
    38  
    39 # disvert : vertical discretisation : string (default='std') : std, ncar 
    40 disvert=ncar 
    41  
    42 # optim_it : mesh optimisation : number of iteration : integer (default=0) 
    43 optim_it=0 
    44  
    45 # initial_state 
    46  
    4753 
    4854# guided_type : string (default=none) : none, ncar 
    4955guided_type=ncar 
    5056 
    51 # etat0 : initial state : string (default=jablonowsky06) : jablonowsky06, academic, ncar 
    52 etat0=ncar 
     57#-------------------- parameters for NCAR test cases ------------------------ 
    5358 
    54 # caldyn : computation type for gcm equation : string (default=gcm) : gcm, adv 
    55 caldyn=adv 
     59# NCAR advection test, initial tracer : string ( default='cos_bell') 
     60# const, slotted_cyl, cos_bell, dbl_cos_bell_q1, dbl_cos_bell_q2, complement, hadley 
     61ncar_adv_shape=complement 
     62 
     63# NCAR advection test, wind field : string (default='deform') : solid, deform, hadley 
     64ncar_adv_wind=deform 
     65 
     66# ncar_dz : model layer thickness in meters: real (default=400)  
     67# used if disvert=ncar 
     68ncar_dz=400 
     69 
     70# ncar_T0 : reference temperature for NCAR test cases : real (default=300)  
     71# also used by disvert if disvert=ncar 
     72ncar_T0=300 
     73 
     74# ncar_p0 : reference pressure for NCAR test cases : real (default=1e5)  
     75# also used by disvert if disvert=ncar 
     76ncar_p0=1e5 
     77 
     78# ncar_disvert_c : exponent for B(eta) : integer (default=1) 
     79# used by disvert if disvert=ncar 
     80ncar_disvert_c=1 

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

    r23 r25  
    133133 
    134134    IF ( iadvtr == iapp_tracvl ) THEN  
     135       PRINT *, 'Advection scheme' 
    135136       bigt = dt*iapp_tracvl 
    136137       DO ind=1,ndomain 

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

    r19 r25  
    2424  TYPE(t_request),POINTER :: req(:)  
    2525  TYPE(t_request),POINTER :: req_u(:)    
    26   PUBLIC :: allocate_caldyn,swap_caldyn,init_williamson,caldyn,write_caldyn,Compute_PV,Compute_enstrophy 
     26  PUBLIC :: allocate_caldyn,swap_caldyn,caldyn,write_caldyn,Compute_PV,Compute_enstrophy 
    2727CONTAINS 
    2828 

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

    r19 r25  
    66  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:) 
    77 
    8 CONTAINS  
     8CONTAINS 
    99!========================================================================= 
    1010 
    … …  
    2121  END SUBROUTINE init_disvert   
    2222 
    23  
    2423  SUBROUTINE disvert(ap,bp,presnivs) 
    2524  USE icosa 
    … …  
    2827  REAL(rstd),INTENT(OUT) :: bp(:) 
    2928  REAL(rstd),INTENT(OUT) :: presnivs(:) 
     29  ! reads from run.def : ncar_dz, ncar_T0, ncar_p0, ncar_disvert_c 
     30  INTEGER :: l,cindx 
     31  REAL(rstd) :: H, eta_top, eta 
     32 
     33  ncar_dz=400 ; CALL getin('ncar_dz',ncar_dz); 
     34 
     35! SELECT CASE(  
     36! pressure profile depends on test case 
     37! coded here for 1-x (transport) 
     38 
     39  ncar_T0=300; CALL getin('ncar_T0',ncar_T0) 
     40  ncar_p0=1e5; CALL getin('ncar_p0',ncar_p0) 
     41  cindx=1 ; CALL getin('ncar_disvert_c',cindx) 
     42 
     43  H = ncar_T0*cpp*kappa/g ! height scale R.T0/g with R=cpp*kappa 
     44  eta_top = exp(-llm*ncar_dz/H) 
     45  
     46  do l = 1,llm+1 
     47     eta = exp(-(l-1)*ncar_dz/H) 
     48     bp(l) = ((eta - eta_top)/(1 - eta_top))**cindx 
     49     ap(l) = ncar_p0 * ( eta - bp(l) ) 
     50  ENDDO 
     51  bp(1)=1. 
     52  ap(1)=0. 
     53  bp(llm+1) = 0 
    3054   
    31   REAL(rstd) :: sig(llm+1) 
    32   REAL(rstd) :: sigtop 
    33   REAL(rstd),PARAMETER:: p0=100000.0  
    34   INTEGER :: l,cindx 
    35   REAL(rstd) :: hdz, ehdz  
     55  DO l = 1, llm 
     56     presnivs(l) = 0.5 *( ap(l)+bp(l)*ncar_p0 + ap(l+1)+bp(l+1)*ncar_p0 ) 
     57  ENDDO 
    3658 
     59  PRINT *, 'Vertical placement of model levels according to DCMIP Appendix E.3' 
     60  PRINT *, 'Parameters : ncar_dz=', ncar_dz, '  ncar_p0=',ncar_p0, '  ncar_disvert_c=',cindx 
     61  PRINT *, 'Isothermal amtosphere with ncar_T0=',ncar_T0  
    3762 
    38          hdz = 400.*g/(300.*287.)  
    39          ehdz = exp(-hdz)  
    40  
    41     do l = 1,llm+1 
    42         sig(l) = ehdz**(l-1)  
    43     end do   
    44          
    45        sigtop = sig(llm+1)  
    46         cindx = 1  
    47  
    48     bp(llm+1) =   0. 
    49     DO l = 1, llm 
    50       bp(l) = (sig(l) - sigtop)/(1 - sigtop) 
    51       bp(l) = bp(l)**cindx 
    52       ap(l) = p0 * ( sig(l) - bp(l) ) 
    53     ENDDO 
    54     bp(1)=1. 
    55     ap(1)=0. 
    56     ap(llm+1) = p0 * ( sig(llm+1) - bp(llm+1) ) 
    57    
    58     DO l = 1, llm 
    59       presnivs(l) = 0.5 *( ap(l)+bp(l)*p0 + ap(l+1)+bp(l+1)*p0 ) 
    60    ENDDO 
    61  
    62   END SUBROUTINE disvert 
     63END SUBROUTINE disvert 
    6364 
    6465END  MODULE disvert_ncar_mod 

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

    r19 r25  
    1919  REAL(rstd), PARAMETER :: rt=radius*0.5 
    2020  REAL(rstd), PARAMETER :: zc=5000.0 
    21   REAL(rstd), PARAMETER :: ps0=100000.0 
    22   REAL(rstd), PARAMETER :: T0=300 
    23   REAL(rstd), PARAMETER :: R_d=287.0    
    2421 
    2522  PUBLIC etat0 
    … …  
    7875  REAL(rstd) :: X2(3),X1(3) 
    7976  INTEGER :: i,j,n,l 
    80   INTEGER :: testcase 
    81  
    82     u = 0.0 ; phis = 0 ; theta_rhodz = 0 ; ps = ps0 
    83          
    84     DO l=1, llm+1 
    85       pr = ap(l) + bp(l)*ps0 
    86       zr(l) = -R_d*T0/g*log(pr/ps0)  
    87     ENDDO     
    88  
    89     DO l=1, llm  
    90       zrl(l) = 0.5*(zr(l) + zr(l+1))  
    91     END DO 
    92  
    93     testcase=5 
    94     CALL getin('ncar_test_case',testcase) 
    95   
    96     SELECT CASE(testcase) 
    97 !--------------------------------------------- SINGLE COSINE BELL 
    98       CASE(0) 
    99         CALL cosine_bell_1(q)  
    100          
    101       CASE(1) 
    102         CALL cosine_bell_2(q)  
    103          
    104       CASE(2) 
    105         CALL cosine_bell_2(q)  
    106         DO l=1,llm 
    107           q(:,l)= 0.9 - 0.8*q(:,l)*q(:,l) 
    108         END DO   
    109          
    110       CASE(3) 
    111         CALL slotted_cylinders(q)        
    112  
    113       CASE(4)  
    114         CALL cosine_bell_2(qxt1) 
    115         DO l = 1,llm 
    116           q(:,l) = 0.9 - 0.8*qxt1(:,l)*qxt1(:,l)   
    117         END DO  
    118         q = q + qxt1  
    119         CALL slotted_cylinders(qxt1) 
    120         q = q + qxt1  
    121         q = 1. - q*0.3   
    122           
    123       CASE(5)  ! hadley like meridional circulation  
    124         CALL hadleyq(q)  
    125          
    126       CASE DEFAULT 
    127         PRINT*,"no such initial profile" 
    128         STOP 
    129          
    130       END SELECT 
    131   
    132   CONTAINS 
     77  CHARACTER(len=255) :: ncar_adv_shape 
     78 
     79  u = 0.0 ; phis = 0 ; theta_rhodz = 0 ; ps = ncar_p0 
     80   
     81  DO l=1, llm+1 
     82     pr = ap(l) + bp(l)*ncar_p0 
     83     zr(l) = -kappa*cpp*ncar_T0/g*log(pr/ncar_p0)  
     84  ENDDO 
     85   
     86  DO l=1, llm  
     87     zrl(l) = 0.5*(zr(l) + zr(l+1))  
     88  END DO 
     89   
     90  ncar_adv_shape='cos_bell' 
     91  CALL getin('ncar_adv_shape',ncar_adv_shape) 
     92   
     93  SELECT CASE(TRIM(ncar_adv_shape)) 
     94     !--------------------------------------------- SINGLE COSINE BELL 
     95  CASE('const') 
     96     q=1 
     97  CASE('cos_bell') 
     98     CALL cosine_bell_1(q)  
     99      
     100  CASE('slotted_cyl') 
     101     CALL slotted_cylinders(q)   
     102      
     103  CASE('dbl_cos_bell_q1') 
     104     CALL cosine_bell_2(q)  
     105      
     106  CASE('dbl_cos_bell_q2') 
     107     CALL cosine_bell_2(q)  
     108     DO l=1,llm 
     109        q(:,l)= 0.9 - 0.8*q(:,l)*q(:,l) 
     110     END DO 
     111      
     112  CASE('complement') 
     113     ! tracer such that, in combination with the other tracer fields  
     114     ! with weight (3/10), the sum is equal to one 
     115     CALL cosine_bell_2(qxt1) 
     116     DO l = 1,llm 
     117        q(:,l) = 0.9 - 0.8*qxt1(:,l)*qxt1(:,l)   
     118     END DO 
     119     q = q + qxt1  
     120     CALL slotted_cylinders(qxt1) 
     121     q = q + qxt1  
     122     q = 1. - q*0.3   
     123      
     124  CASE('hadley')  ! hadley like meridional circulation  
     125     CALL hadleyq(q)  
     126      
     127  CASE DEFAULT 
     128     PRINT *, 'Bad selector for variable ncar_adv_shape : <', TRIM(ncar_adv_shape),   & 
     129          '> options are <const>, <slotted_cyl>, <cos_bell>, <dbl_cos_bell_q1>', & 
     130          '<dbl_cos_bell_q2>, <complement>, <hadley>' 
     131     STOP 
     132      
     133  END SELECT 
     134       
     135CONTAINS 
    133136 
    134137!====================================================================== 

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

    r19 r25  
    1919       
    2020      CASE ('ncar') 
    21         CALL init_guided_ncar(dt) 
     21        CALL init_guided_ncar 
    2222         
    2323      CASE DEFAULT 
    … …  
    3030 
    3131   
    32   SUBROUTINE guided(it, f_ps, f_theta_rhodz, f_u, f_q) 
     32  SUBROUTINE guided(tt, f_ps, f_theta_rhodz, f_u, f_q) 
    3333  USE icosa 
    3434  USE guided_ncar_mod, ONLY : guided_ncar => guided 
    3535  IMPLICIT NONE 
    36     INTEGER, INTENT(IN)   :: it 
     36    REAL(rstd), INTENT(IN):: tt 
    3737    TYPE(t_field),POINTER :: f_ps(:) 
    3838    TYPE(t_field),POINTER :: f_phis(:) 
    … …  
    4545       
    4646      CASE ('ncar') 
    47         CALL guided_ncar(it, f_ps, f_theta_rhodz, f_u, f_q) 
     47        CALL guided_ncar(tt, f_ps, f_theta_rhodz, f_u, f_q) 
    4848    END SELECT  
    4949   

  • 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 

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

    r19 r25  
    1515  USE transfert_mod 
    1616   
     17  ! Variables defined by run.def 
     18 
     19  REAL(rstd) :: ncar_dz, ncar_p0, ncar_T0 ! read from run.def by disvert 
     20 
    1721END MODULE icosa 

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

    r21 r25  
    9090  CALL init_caldyn(dt) 
    9191  CALL init_guided(dt) 
    92 !  CALL init_advect_tracer(dt) 
     92  CALL init_advect_tracer(dt) 
    9393   
    9494  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    … …  
    9797    PRINT *,"It No :",It 
    9898 
    99     CALL guided(it,f_ps,f_theta_rhodz,f_u,f_q) 
     99    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) 
    100100    CALL caldyn(it,f_phis,f_ps,f_theta_rhodz,f_u, f_dps, f_dtheta_rhodz, f_du) 
    101 !    CALL advect_tracer(f_ps,f_u,f_q) 
     101    CALL advect_tracer(f_ps,f_u,f_q) 
    102102     
    103103    SELECT CASE (TRIM(scheme)) 
Note: See TracChangeset for help on using the changeset viewer.