Changeset 25 for codes/icosagcm/trunk
- Timestamp:
- 07/26/12 04:05:44 (12 years ago)
- Location:
- codes/icosagcm/trunk
- Files:
-
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/run_adv.def
r18 r25 1 #------------------------------- Mesh --------------------------------- 2 1 3 # Number of subdivision on a main triangle (nbp) : integer (default=40) 2 4 nbp=20 3 5 4 # Number of vertical layer (llm) : integer (default=19) 6 # optim_it : mesh optimisation : number of iteration : integer (default=0) 7 optim_it=0 8 9 # Number of vertical layers (llm) : integer (default=19) 5 10 llm=19 11 12 # disvert : vertical discretisation : string (default='std') : std, ncar 13 disvert=ncar 14 15 #---------------------------------- Misc -------------------------------- 6 16 7 17 # number of tracer (nqtot) : integer (default 1) … … 11 21 dt = 180. 12 22 13 # scheme type : string ( default='adam_bashforth') euler, leapfrog, leapfrog_matsuno, adam_bashforth) 23 # number of timestep (default 100) 24 itaumax = 38400 25 26 # output field period : integer (default none) 27 write_period=3600 28 29 # etat0 : initial state : string (default=jablonowsky06) : 30 # jablonowsky06, academic, ncar 31 etat0=ncar 32 33 # caldyn : computation type for gcm equation : string (default=gcm) : gcm, adv 34 caldyn=adv 35 36 #------------------------------ Dynamics -------------------------------- 37 38 # scheme type : string ( default='adam_bashforth') 39 # euler, leapfrog, leapfrog_matsuno, adam_bashforth) 14 40 scheme = euler 15 41 16 42 # matsuno period : integer ( default=5) 17 43 matsuno_period = 5 18 19 # number of timestep (default 100)20 itaumax = 3840021 44 22 45 # dissipation time graddiv : real (default=5000) … … 25 48 26 49 # dissipation time nxgradrot (default=5000) 27 tetagrot = 1800 28 tau_gradrot = 1800 29 50 tetagrot=1800 51 tau_gradrot=1800 30 52 tau_divgrad=1800 31 32 # output field period : integer (default none)33 write_period=360034 35 36 #NCAR testcase : integer ( default=5)37 ncar_test_case=138 39 # disvert : vertical discretisation : string (default='std') : std, ncar40 disvert=ncar41 42 # optim_it : mesh optimisation : number of iteration : integer (default=0)43 optim_it=044 45 # initial_state46 47 53 48 54 # guided_type : string (default=none) : none, ncar 49 55 guided_type=ncar 50 56 51 # etat0 : initial state : string (default=jablonowsky06) : jablonowsky06, academic, ncar 52 etat0=ncar 57 #-------------------- parameters for NCAR test cases ------------------------ 53 58 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 61 ncar_adv_shape=complement 62 63 # NCAR advection test, wind field : string (default='deform') : solid, deform, hadley 64 ncar_adv_wind=deform 65 66 # ncar_dz : model layer thickness in meters: real (default=400) 67 # used if disvert=ncar 68 ncar_dz=400 69 70 # ncar_T0 : reference temperature for NCAR test cases : real (default=300) 71 # also used by disvert if disvert=ncar 72 ncar_T0=300 73 74 # ncar_p0 : reference pressure for NCAR test cases : real (default=1e5) 75 # also used by disvert if disvert=ncar 76 ncar_p0=1e5 77 78 # ncar_disvert_c : exponent for B(eta) : integer (default=1) 79 # used by disvert if disvert=ncar 80 ncar_disvert_c=1 -
codes/icosagcm/trunk/src/advect_tracer.f90
r23 r25 133 133 134 134 IF ( iadvtr == iapp_tracvl ) THEN 135 PRINT *, 'Advection scheme' 135 136 bigt = dt*iapp_tracvl 136 137 DO ind=1,ndomain -
codes/icosagcm/trunk/src/caldyn_sw.f90
r19 r25 24 24 TYPE(t_request),POINTER :: req(:) 25 25 TYPE(t_request),POINTER :: req_u(:) 26 PUBLIC :: allocate_caldyn,swap_caldyn, init_williamson,caldyn,write_caldyn,Compute_PV,Compute_enstrophy26 PUBLIC :: allocate_caldyn,swap_caldyn,caldyn,write_caldyn,Compute_PV,Compute_enstrophy 27 27 CONTAINS 28 28 -
codes/icosagcm/trunk/src/disvert_ncar.f90
r19 r25 6 6 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:) 7 7 8 CONTAINS 8 CONTAINS 9 9 !========================================================================= 10 10 … … 21 21 END SUBROUTINE init_disvert 22 22 23 24 23 SUBROUTINE disvert(ap,bp,presnivs) 25 24 USE icosa … … 28 27 REAL(rstd),INTENT(OUT) :: bp(:) 29 28 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 30 54 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 36 58 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 37 62 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 63 END SUBROUTINE disvert 63 64 64 65 END MODULE disvert_ncar_mod -
codes/icosagcm/trunk/src/etat0_ncar.f90
r19 r25 19 19 REAL(rstd), PARAMETER :: rt=radius*0.5 20 20 REAL(rstd), PARAMETER :: zc=5000.0 21 REAL(rstd), PARAMETER :: ps0=100000.022 REAL(rstd), PARAMETER :: T0=30023 REAL(rstd), PARAMETER :: R_d=287.024 21 25 22 PUBLIC etat0 … … 78 75 REAL(rstd) :: X2(3),X1(3) 79 76 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 135 CONTAINS 133 136 134 137 !====================================================================== -
codes/icosagcm/trunk/src/guided_mod.f90
r19 r25 19 19 20 20 CASE ('ncar') 21 CALL init_guided_ncar (dt)21 CALL init_guided_ncar 22 22 23 23 CASE DEFAULT … … 30 30 31 31 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) 33 33 USE icosa 34 34 USE guided_ncar_mod, ONLY : guided_ncar => guided 35 35 IMPLICIT NONE 36 INTEGER, INTENT(IN) :: it36 REAL(rstd), INTENT(IN):: tt 37 37 TYPE(t_field),POINTER :: f_ps(:) 38 38 TYPE(t_field),POINTER :: f_phis(:) … … 45 45 46 46 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) 48 48 END SELECT 49 49 -
codes/icosagcm/trunk/src/guided_ncar_mod.f90
r19 r25 3 3 PRIVATE 4 4 5 REAL(rstd),SAVE :: dt 6 INTEGER,SAVE :: test 5 INTEGER,SAVE :: case_wind 7 6 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 8 13 PUBLIC init_guided, guided 9 14 10 15 CONTAINS 11 16 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 20 32 END SUBROUTINE init_guided 21 22 33 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) 24 35 USE icosa 25 36 IMPLICIT NONE 26 INTEGER, INTENT(IN) :: it37 REAL(rstd), INTENT(IN):: tt 27 38 TYPE(t_field),POINTER :: f_ps(:) 28 39 TYPE(t_field),POINTER :: f_phis(:) … … 38 49 CALL swap_geometry(ind) 39 50 ue = f_u(ind) 40 CALL wind_profile( it,ue)51 CALL wind_profile(tt,ue) 41 52 END DO 42 53 … … 44 55 45 56 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 50 62 REAL(rstd),INTENT(OUT) :: ue(iim*3*jjm,llm) 51 63 REAL(rstd) :: lon, lat 52 64 REAL(rstd) :: nx(3),n_norm,Velocity(3,llm) 53 REAL(rstd), PARAMETER :: lon0=3*pi/2,lat0=0.054 65 REAL(rstd) :: rr1,rr2,bb,cc,aa,hmx 55 66 REAL(rstd) :: v1(3),v2(3),ny(3) 56 67 INTEGER :: i,j,n,l 57 REAL(rstd) :: zrl(llm)68 REAL(rstd) :: pitbytau,kk, pr, zr, u0, u1, v0 58 69 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 60 73 !--------------------------------------------------------- 61 74 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 74 92 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-1079 n_norm=sqrt(sum(nx(:)**2))80 IF (n_norm>1e-30) THEN81 nx=-nx/n_norm*ne(n,lup)82 ue(n+u_lup,l)=sum(nx(:)*velocity(:,l))83 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-1089 n_norm=sqrt(sum(nx(:)**2))90 IF (n_norm>1e-30) THEN91 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,j94 ENDIF95 ENDDO96 ENDDO93 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 97 115 END DO 98 116 99 117 CONTAINS 100 118 … … 104 122 INTEGER,INTENT(IN)::l 105 123 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) 108 125 REAL(rstd) :: lon,lat 109 REAL(rstd),PARAMETER::alpha=0.0110 126 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 125 128 CALL xyz2lonlat(x/radius,lon,lat) 126 129 e_lat(1) = -cos(lon)*sin(lat) … … 133 136 134 137 u = 0.0 ; v = 0.0 135 u0 = 2*pi*radius/tau136 138 137 SELECT CASE( test)138 CASE(0)139 140 141 CASE(1)142 lon = lon - 2*pitbytau143 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 150 151 END SELECT152 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 153 155 Velocity(:)=(u*e_lon(:)+v*e_lat(:)+1e-50) 154 156 155 157 END SUBROUTINE compute_velocity 156 158 157 SUBROUTINE compute_zrl(zrl)158 USE disvert_mod159 IMPLICIT NONE160 REAL(rstd),INTENT(OUT) :: zrl(llm)161 INTEGER :: l162 REAL(rstd) :: zr(llm+1)163 REAl(rstd) :: pr164 REAL(rstd), PARAMETER :: T0=300165 REAL(rstd), PARAMETER :: R_d=287.0166 REAL(rstd), PARAMETER :: ps0=100000.0167 168 169 DO l = 1, llm+1170 pr = ap(l) + bp(l)*ps0171 zr(l) = -R_d*T0/g*log(pr/ps0)172 ENDDO173 174 DO l = 1, llm175 zrl(l) = 0.5*(zr(l) + zr(l+1))176 END DO177 178 END SUBROUTINE compute_zrl179 180 159 END SUBROUTINE wind_profile 181 160 -
codes/icosagcm/trunk/src/icosa_mod.f90
r19 r25 15 15 USE transfert_mod 16 16 17 ! Variables defined by run.def 18 19 REAL(rstd) :: ncar_dz, ncar_p0, ncar_T0 ! read from run.def by disvert 20 17 21 END MODULE icosa -
codes/icosagcm/trunk/src/timeloop_gcm.f90
r21 r25 90 90 CALL init_caldyn(dt) 91 91 CALL init_guided(dt) 92 !CALL init_advect_tracer(dt)92 CALL init_advect_tracer(dt) 93 93 94 94 CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) … … 97 97 PRINT *,"It No :",It 98 98 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) 100 100 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) 102 102 103 103 SELECT CASE (TRIM(scheme))
Note: See TracChangeset
for help on using the changeset viewer.