Changeset 129 for codes/icosagcm/trunk
- Timestamp:
- 02/04/13 15:15:35 (12 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/caldyn.f90
r110 r129 30 30 END SUBROUTINE init_caldyn 31 31 32 SUBROUTINE caldyn( it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du)32 SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 33 33 USE icosa 34 34 USE caldyn_gcm_mod, ONLY : caldyn_gcm=>caldyn 35 35 USE caldyn_adv_mod, ONLY : caldyn_adv=>caldyn 36 36 IMPLICIT NONE 37 INTEGER,INTENT(IN) :: it37 LOGICAL,INTENT(IN) :: write_out 38 38 TYPE(t_field),POINTER :: f_phis(:) 39 39 TYPE(t_field),POINTER :: f_ps(:) … … 47 47 SELECT CASE (TRIM(caldyn_type)) 48 48 CASE('gcm') 49 CALL caldyn_gcm( it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du)49 CALL caldyn_gcm(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 50 50 CASE('adv') 51 CALL caldyn_adv( it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du)51 CALL caldyn_adv(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 52 52 END SELECT 53 53 -
codes/icosagcm/trunk/src/caldyn_adv.f90
r110 r129 79 79 80 80 81 SUBROUTINE caldyn( it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du)81 SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 82 82 USE icosa 83 83 USE vorticity_mod … … 85 85 USE theta2theta_rhodz_mod 86 86 IMPLICIT NONE 87 INTEGER,INTENT(IN) :: it87 LOGICAL,INTENT(IN) :: write_out 88 88 TYPE(t_field),POINTER :: f_phis(:) 89 89 TYPE(t_field),POINTER :: f_ps(:) … … 133 133 134 134 135 IF ( mod(it,itau_out)==0) THEN135 IF (write_out) THEN 136 136 CALL writefield("ps",f_ps) 137 137 ENDIF -
codes/icosagcm/trunk/src/caldyn_gcm.f90
r128 r129 80 80 END SUBROUTINE allocate_caldyn 81 81 82 SUBROUTINE caldyn( it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du)82 SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 83 83 USE icosa 84 84 USE vorticity_mod … … 86 86 USE theta2theta_rhodz_mod 87 87 IMPLICIT NONE 88 INTEGER,INTENT(IN) :: it88 LOGICAL,INTENT(IN) :: write_out 89 89 TYPE(t_field),POINTER :: f_phis(:) 90 90 TYPE(t_field),POINTER :: f_ps(:) … … 168 168 END SELECT 169 169 170 IF ( mod(it,itau_out)==0) THEN170 IF (write_out) THEN 171 171 PRINT *,'CALL write_output_fields' 172 172 CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, & -
codes/icosagcm/trunk/src/dissip_gcm.f90
r109 r129 250 250 ENDDO 251 251 PRINT *,"gradiv : dumax",dumax 252 252 PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., & 253 'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000 254 PRINT *, 'Max. time step assuming c=340 m/s and Courant number=2.8 :', & 255 2.8/340.*dumax**(-.5/nitergdiv) 256 253 257 cgraddiv=dumax**(-1./nitergdiv) 254 258 PRINT *,"cgraddiv : ",cgraddiv … … 469 473 n=(j-1)*iim+i 470 474 471 due(n+u_right,l)=due(n+u_right,l)-0.5*(due_diss1(n+u_right,l)/tau_graddiv(l)+due_diss2(n+u_right,l)/tau_gradrot(l)) 472 due(n+u_lup,l)=due(n+u_lup,l) -0.5*(due_diss1(n+u_lup,l) /tau_graddiv(l)+due_diss2(n+u_lup,l) /tau_gradrot(l)) 473 due(n+u_ldown,l)=due(n+u_ldown,l)-0.5*(due_diss1(n+u_ldown,l)/tau_graddiv(l)+due_diss2(n+u_ldown,l)/tau_gradrot(l)) 474 475 dtheta_rhodz(n,l)=dtheta_rhodz(n,l)-0.5*dtheta_rhodz_diss(n,l)/tau_divgrad(l) 476 475 due(n+u_right,l) = -0.5*( due_diss1(n+u_right,l)/tau_graddiv(l) + due_diss2(n+u_right,l)/tau_gradrot(l) ) 476 due(n+u_lup,l) = -0.5*( due_diss1(n+u_lup,l) /tau_graddiv(l) + due_diss2(n+u_lup,l) /tau_gradrot(l) ) 477 due(n+u_ldown,l) = -0.5*( due_diss1(n+u_ldown,l)/tau_graddiv(l) + due_diss2(n+u_ldown,l)/tau_gradrot(l) ) 478 479 dtheta_rhodz(n,l) = -0.5*dtheta_rhodz_diss(n,l)/tau_divgrad(l) 477 480 ENDDO 478 481 ENDDO -
codes/icosagcm/trunk/src/timeloop_gcm.f90
r124 r129 14 14 IMPLICIT NONE 15 15 TYPE(t_field),POINTER :: f_phis(:) 16 TYPE(t_field),POINTER :: f_theta(:)16 ! TYPE(t_field),POINTER :: f_theta(:) 17 17 TYPE(t_field),POINTER :: f_q(:) 18 18 TYPE(t_field),POINTER :: f_dtheta(:) … … 35 35 ! REAL(rstd) :: dt, run_length 36 36 INTEGER :: ind 37 INTEGER :: it,i,j,n 37 INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period 38 38 CHARACTER(len=255) :: scheme 39 INTEGER :: matsuno_period40 39 ! INTEGER :: itaumax 41 40 ! REAL(rstd) ::write_period … … 54 53 ! dt=dt/scale_factor 55 54 56 scheme='adam_bashforth' 55 ! Trends 56 CALL allocate_field(f_dps,field_t,type_real) 57 CALL allocate_field(f_du,field_u,type_real,llm) 58 CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm) 59 ! Model state at current time step (RK/MLF/Euler) 60 CALL allocate_field(f_phis,field_t,type_real) 61 CALL allocate_field(f_ps,field_t,type_real) 62 CALL allocate_field(f_u,field_u,type_real,llm) 63 CALL allocate_field(f_theta_rhodz,field_t,type_real,llm) 64 ! Model state at previous time step (RK/MLF) 65 CALL allocate_field(f_psm1,field_t,type_real) 66 CALL allocate_field(f_um1,field_u,type_real,llm) 67 CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm) 68 ! Tracers 69 CALL allocate_field(f_q,field_t,type_real,llm,nqtot) 70 71 scheme='runge_kutta' 57 72 CALL getin('scheme',scheme) 58 59 matsuno_period=5 60 CALL getin('matsuno_period',matsuno_period) 61 IF (TRIM(scheme)=='leapfrog') matsuno_period=itaumax+1 62 73 74 SELECT CASE (TRIM(scheme)) 75 CASE('euler') 76 nb_stage=1 77 CASE ('runge_kutta') 78 nb_stage=4 79 CASE ('leapfrog_matsuno') 80 matsuno_period=5 81 CALL getin('matsuno_period',matsuno_period) 82 nb_stage=matsuno_period+1 83 ! Model state 2 time steps ago (MLF) 84 CALL allocate_field(f_psm2,field_t,type_real) 85 CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm) 86 CALL allocate_field(f_um2,field_u,type_real,llm) 87 CASE default 88 PRINT*,'Bad selector for variable scheme : <', TRIM(scheme), & 89 ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>' 90 STOP 91 END SELECT 92 63 93 ! write_period=0 64 94 ! CALL getin('write_period',write_period) … … 67 97 ! PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out 68 98 69 CALL allocate_field(f_phis,field_t,type_real) 70 71 CALL allocate_field(f_ps,field_t,type_real) 72 CALL allocate_field(f_psm1,field_t,type_real) 73 CALL allocate_field(f_psm2,field_t,type_real) 74 CALL allocate_field(f_dps,field_t,type_real) 75 CALL allocate_field(f_dpsm1,field_t,type_real) 76 CALL allocate_field(f_dpsm2,field_t,type_real) 77 78 CALL allocate_field(f_u,field_u,type_real,llm) 79 CALL allocate_field(f_um1,field_u,type_real,llm) 80 CALL allocate_field(f_um2,field_u,type_real,llm) 81 CALL allocate_field(f_du,field_u,type_real,llm) 82 CALL allocate_field(f_dum1,field_u,type_real,llm) 83 CALL allocate_field(f_dum2,field_u,type_real,llm) 84 85 CALL allocate_field(f_theta,field_t,type_real,llm) 86 CALL allocate_field(f_dtheta,field_t,type_real,llm) 87 88 CALL allocate_field(f_q,field_t,type_real,llm,nqtot) 89 90 CALL allocate_field(f_theta_rhodz,field_t,type_real,llm) 91 CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm) 92 CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm) 93 CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm) 94 CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm) 95 CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm) 99 ! Trends at previous time steps needed only by Adams-Bashforth 100 ! CALL allocate_field(f_dpsm1,field_t,type_real) 101 ! CALL allocate_field(f_dpsm2,field_t,type_real) 102 ! CALL allocate_field(f_dum1,field_u,type_real,llm) 103 ! CALL allocate_field(f_dum2,field_u,type_real,llm) 104 ! CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm) 105 ! CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm) 106 107 ! CALL allocate_field(f_theta,field_t,type_real,llm) 108 ! CALL allocate_field(f_dtheta,field_t,type_real,llm) 96 109 97 110 CALL init_dissip … … 113 126 ENDIF 114 127 115 CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) 116 CALL caldyn(it,f_phis,f_ps,f_theta_rhodz,f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 128 ! CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) 129 130 DO stage=1,nb_stage 131 CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), & 132 f_phis,f_ps,f_theta_rhodz,f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 133 SELECT CASE (TRIM(scheme)) 134 CASE('euler') 135 CALL euler_scheme 136 137 CASE ('runge_kutta') 138 CALL rk_scheme(stage) 139 140 CASE ('leapfrog_matsuno') 141 CALL leapfrog_matsuno_scheme(stage) 142 143 ! CASE ('leapfrog') 144 ! CALL leapfrog_scheme 145 ! 146 ! CASE ('adam_bashforth') 147 ! CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz) 148 ! CALL adam_bashforth_scheme 149 CASE default 150 ! PRINT*,'Bad selector for variable scheme : <', TRIM(scheme), & 151 ! ' > options are <euler>, <runge_kutta>, <leapfrog>, <leapfrog_matsuno>, <adam_bashforth>' 152 PRINT*,'Bad selector for variable scheme : <', TRIM(scheme), & 153 ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>' 154 STOP 155 156 END SELECT 157 END DO 117 158 ! CALL advect_tracer(f_ps,f_u,f_q) 118 119 SELECT CASE (TRIM(scheme)) 120 CASE('euler') 121 CALL euler_scheme 122 123 CASE ('leapfrog') 124 CALL leapfrog_scheme 125 126 CASE ('runge_kutta') 127 CALL rk_scheme 128 129 CASE ('leapfrog_matsuno') 130 CALL leapfrog_matsuno_scheme 131 132 CASE ('adam_bashforth') 133 CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz) 134 CALL adam_bashforth_scheme 135 136 CASE default 137 PRINT*,'Bad selector for variable scheme : <', TRIM(scheme), & 138 ' > options are <euler>, <runge_kutta>, <leapfrog>, <leapfrog_matsuno>, <adam_bashforth>' 139 STOP 140 141 END SELECT 142 159 ! CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz) 143 160 ! CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q) 144 145 ENDDO 146 161 ENDDO 162 147 163 CONTAINS 148 164 … … 161 177 END SUBROUTINE Euler_scheme 162 178 163 SUBROUTINE RK_scheme 179 SUBROUTINE RK_scheme(stage) 164 180 IMPLICIT NONE 165 181 INTEGER :: ind, stage 166 REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ 1., 4./3., 2., 4. /) 182 ! REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ 1., 4./3., 2., 4. /) 183 REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /) 167 184 REAL(rstd) :: tau 168 185 169 stage = MOD(it,4)+1170 186 tau = dt*coef(stage) 171 187 172 188 DO ind=1,ndomain 173 189 ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) 174 190 psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind) 175 191 dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind) 176 177 IF (stage==1) THEN 192 193 IF (stage==1) THEN ! first stage : save model state in XXm1 178 194 psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) 179 195 END IF 196 ! updates are of the form x1 := x0 + tau*f(x1) 180 197 ps(:)=psm1(:)+tau*dps(:) 181 198 u(:,:)=um1(:,:)+tau*du(:,:) 182 199 theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:) 183 200 END DO 201 184 202 END SUBROUTINE RK_scheme 185 203 … … 215 233 END SUBROUTINE leapfrog_scheme 216 234 217 SUBROUTINE leapfrog_matsuno_scheme 235 SUBROUTINE leapfrog_matsuno_scheme(stage) 218 236 IMPLICIT NONE 219 INTEGER :: ind 220 237 INTEGER :: ind, stage 238 REAL :: tau 239 tau = dt/nb_stage 221 240 DO ind=1,ndomain 222 241 ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) … … 226 245 227 246 228 IF (MOD(it,matsuno_period+1)==0) THEN 247 ! IF (MOD(it,matsuno_period+1)==0) THEN 248 IF (stage==1) THEN 229 249 psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) 230 ps(:)=psm1(:)+dt*dps(:) 231 u(:,:)=um1(:,:)+dt*du(:,:) 232 theta_rhodz(:,:)=theta_rhodzm1(:,:)+dt*dtheta_rhodz(:,:) 233 234 ELSE IF (MOD(it,matsuno_period+1)==1) THEN 235 236 ps(:)=psm1(:)+dt*dps(:) 237 u(:,:)=um1(:,:)+dt*du(:,:) 238 theta_rhodz(:,:)=theta_rhodzm1(:,:)+dt*dtheta_rhodz(:,:) 250 ps(:)=psm1(:)+tau*dps(:) 251 u(:,:)=um1(:,:)+tau*du(:,:) 252 theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:) 253 254 ! ELSE IF (MOD(it,matsuno_period+1)==1) THEN 255 ELSE IF (stage==2) THEN 256 257 ps(:)=psm1(:)+tau*dps(:) 258 u(:,:)=um1(:,:)+tau*du(:,:) 259 theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:) 239 260 240 261 psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) … … 243 264 ELSE 244 265 245 ps(:)=psm2(:)+2* dt*dps(:)246 u(:,:)=um2(:,:)+2* dt*du(:,:)247 theta_rhodz(:,:)=theta_rhodzm2(:,:)+2* dt*dtheta_rhodz(:,:)266 ps(:)=psm2(:)+2*tau*dps(:) 267 u(:,:)=um2(:,:)+2*tau*du(:,:) 268 theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:) 248 269 249 270 psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)
Note: See TracChangeset
for help on using the changeset viewer.