[4] | 1 | !> \file step.F90 |
---|
| 2 | !! TOOOOOOOOO DOOOOOOOOOOOOO |
---|
| 3 | !! Step_grisli contient |
---|
| 4 | !< |
---|
| 5 | |
---|
| 6 | |
---|
| 7 | subroutine step_grisli() |
---|
| 8 | |
---|
| 9 | use module3d_phy |
---|
| 10 | use module_choix ! module de choix du type de run |
---|
| 11 | ! module_choix donne acces a tous les modules |
---|
| 12 | ! de declaration des packages |
---|
| 13 | use interface_icetempmod |
---|
| 14 | use sorties_ncdf_grisli |
---|
| 15 | use flottab_mod |
---|
| 16 | use diagno_mod |
---|
| 17 | use resolmeca_SIA_L1 |
---|
| 18 | use track_debug |
---|
| 19 | |
---|
| 20 | use geom_type |
---|
| 21 | use temperature_type |
---|
| 22 | use ice_flow_type |
---|
| 23 | use mask_flgz_type |
---|
| 24 | use deformation_type |
---|
| 25 | use autre_pr_temp_type |
---|
| 26 | use step_temp_type |
---|
| 27 | |
---|
| 28 | |
---|
| 29 | implicit none |
---|
| 30 | |
---|
| 31 | integer :: m |
---|
| 32 | logical :: shelf_vitbil = .true. ! si vrai on prend les vitesses de bilan pour le shelf |
---|
| 33 | |
---|
| 34 | !!$if (time.lt.timemax-100) then |
---|
| 35 | !!$ print*,'le modele remonte le temps de',timemax,'a',time |
---|
| 36 | !!$ stop |
---|
| 37 | !!$endif |
---|
| 38 | timemax=time |
---|
| 39 | |
---|
| 40 | if (itracebug.eq.1) write(num_tracebug,*) 'entree dans step_grisli', & |
---|
| 41 | ' nt=',nt,'iter_beta=',iter_beta |
---|
| 42 | !-------------------------------------------------------------------- |
---|
| 43 | ! lecture dans un petit fichier au cas ou les parametres du run |
---|
| 44 | ! seraient changes |
---|
| 45 | !-------------------------------------------------------------------- |
---|
| 46 | |
---|
| 47 | call limit_file(nt,real(time),dt,tend,dtsortie,dtcpt,testdiag,dtt,runname) |
---|
| 48 | |
---|
| 49 | if (time.lt.tgrounded) then |
---|
| 50 | marine=.false. |
---|
| 51 | else |
---|
| 52 | marine=.true. |
---|
| 53 | endif |
---|
| 54 | |
---|
| 55 | !if ((iter_beta.eq.0).or.(nt.le.1)) then |
---|
| 56 | call masque() |
---|
| 57 | call flottab() |
---|
| 58 | if (iter_beta.gt.0) call init_dragging |
---|
| 59 | !end if |
---|
| 60 | |
---|
| 61 | |
---|
| 62 | !-----------------------------------------------------------------------! |
---|
| 63 | ! debut bloc appels tous les dtt ! |
---|
| 64 | !------------- -------------------------------------------------------! |
---|
| 65 | |
---|
| 66 | |
---|
| 67 | isync_1: if (ISYNCHRO.eq.1) then !debut bloc appels tous les dtt |
---|
| 68 | |
---|
| 69 | |
---|
| 70 | |
---|
| 71 | spinup_1_veloc: if (ispinup.ne.2.or.nt.lt.5) then |
---|
| 72 | |
---|
| 73 | !--------------------------------------------------------------------------------- |
---|
| 74 | ! spinup_1 : calcul des vitesses dynamiques : |
---|
| 75 | !------------ |
---|
| 76 | ! - est fait dans tous les cas pour les premieres boucles temporelles |
---|
| 77 | ! - est fait la plupart des types de spinup sauf ispinup=2 |
---|
| 78 | ! - ispinup=1 -> modif: maintenant on fait le calcul de veloc |
---|
| 79 | !--------------------------------------------------------------------------------- |
---|
| 80 | |
---|
| 81 | if (itracebug.eq.1) call tracebug(' Dans spinup_1_veloc : calculer velocities') |
---|
| 82 | |
---|
| 83 | call forclim() ! |
---|
| 84 | call ablation |
---|
| 85 | |
---|
| 86 | call Neffect() |
---|
| 87 | call velocities() |
---|
| 88 | endif spinup_1_veloc |
---|
| 89 | |
---|
| 90 | if (nt.gt.2) then |
---|
| 91 | ! Hassine ca va etre modifier bien sur apres: l appel sera beaucoup moin long que ca |
---|
| 92 | |
---|
| 93 | call Allocate_geom(geom_g,nx,ny) |
---|
| 94 | call Init_geom(Geom_g) |
---|
| 95 | |
---|
| 96 | call Allocate_temperature(temperature_g,nx,ny,nz,nz+nzm) |
---|
| 97 | call Init_Temperature(Temperature_g) |
---|
| 98 | |
---|
| 99 | call Allocate_Ice_flow(Ice_flow_g,nx,ny,nz) |
---|
| 100 | call Init_Ice_flow(Ice_flow_g) |
---|
| 101 | |
---|
| 102 | call Allocate_Mask_flgz(Mask_flgz_g,nx,ny) |
---|
| 103 | call Init_Mask_flgz(Mask_flgz_g) |
---|
| 104 | |
---|
| 105 | call Allocate_Deformation(Deformation_g,nx,ny,nz,n1poly,n2poly) |
---|
| 106 | call Init_Deformation_h(Deformation_g) |
---|
| 107 | |
---|
| 108 | call Allocate_autre_pr_temp(Autre_pr_temp_g,nx,ny) |
---|
| 109 | call Init_autre_pr_temp(Autre_pr_temp_g) |
---|
| 110 | |
---|
| 111 | call Init_step_temp(step_temp_g) |
---|
| 112 | |
---|
| 113 | call icetemp(geom_g,temperature_g, & |
---|
| 114 | Ice_flow_g,Mask_flgz_g,& |
---|
| 115 | deformation_g,Autre_pr_temp_g,& |
---|
| 116 | step_temp_g,num_tracebug) |
---|
| 117 | |
---|
| 118 | T= Temperature_g%Temperature |
---|
| 119 | Tpmp= Temperature_g%Temp_melting |
---|
| 120 | Tbdot= Autre_pr_temp_g%Tbdot |
---|
| 121 | Bmelt= Autre_pr_temp_g%BMELT |
---|
| 122 | Ibase=Autre_pr_temp_g%IBASE |
---|
| 123 | Phid=Autre_pr_temp_g%Phid |
---|
| 124 | |
---|
| 125 | ! write(6,*) time,'call icetemp' |
---|
| 126 | ! subroutines pour le calcul de la fusion basale |
---|
| 127 | call bmeltshelf |
---|
| 128 | call bmelt_grounded |
---|
| 129 | |
---|
| 130 | |
---|
| 131 | endif |
---|
| 132 | |
---|
| 133 | spinup_2_flowlaw: if ((ispinup.ne.2).or.(nt.lt.5)) then |
---|
| 134 | |
---|
| 135 | !--------------------------------------------------------------------------------- |
---|
| 136 | ! spinup_2_flowlaw |
---|
| 137 | !---------------- |
---|
| 138 | ! si ispinup = 2 , on ne fait pas le calcul de flowlaw, ni de diffusiv |
---|
| 139 | ! puisque les vitesses de bilan sont imposees (initialisee nt<5) |
---|
| 140 | ! |
---|
| 141 | ! Flowlaw et diffusiv : servent soit directement (ispinup = 0,1) soit pour le |
---|
| 142 | ! rescaling (ispinup = 3) |
---|
| 143 | !--------------------------------------------------------------------------------- |
---|
| 144 | |
---|
| 145 | if (itracebug.eq.1) call tracebug(' Dans spinup_2_flowlaw') |
---|
| 146 | |
---|
| 147 | if (ispinup.eq.3) then |
---|
| 148 | call diffusiv ! pour vérifier les champs |
---|
| 149 | debug_3D(:,:,71)=ux(:,:,1) |
---|
| 150 | debug_3D(:,:,72)=uy(:,:,1) |
---|
| 151 | debug_3D(:,:,13)=ux(:,:,nz) |
---|
| 152 | debug_3D(:,:,14)=uy(:,:,nz) |
---|
| 153 | end if |
---|
| 154 | |
---|
| 155 | call flow_general |
---|
| 156 | do iglen=n1poly,n2poly |
---|
| 157 | call flowlaw(iglen) |
---|
| 158 | end do |
---|
| 159 | |
---|
| 160 | |
---|
| 161 | ! debug_3D(:,:,65)=s2a_mx(:,:,1,1) |
---|
| 162 | ! debug_3D(:,:,66)=s2a_mx(:,:,1,2) |
---|
| 163 | ! debug_3D(:,:,67)=s2a_my(:,:,1,1) |
---|
| 164 | ! debug_3D(:,:,68)=s2a_my(:,:,1,2) |
---|
| 165 | ! debug_3D(:,:,69)=ddx(:,:,1) |
---|
| 166 | ! debug_3D(:,:,69)=ddx(:,:,2) |
---|
| 167 | |
---|
| 168 | if (ispinup.eq.3) then |
---|
| 169 | call calc_coef_vitbil ! recalibre sa_mx et s2a_mx et ddbx |
---|
| 170 | ! call limit_coef_vitbil ! garde la loi de deformation (debug) |
---|
| 171 | end if |
---|
| 172 | |
---|
| 173 | end if spinup_2_flowlaw |
---|
| 174 | |
---|
| 175 | endif isync_1 |
---|
| 176 | |
---|
| 177 | |
---|
| 178 | spinup_3_bed: if (ispinup.eq.0.or.nt.lt.5) then |
---|
| 179 | !--------------------------------------------------------------------------------- |
---|
| 180 | ! spinup_3_bed |
---|
| 181 | !---------------- |
---|
| 182 | ! Ne passe dans ce bloc que quand on veut l'isostasie, donc si variations epaisseur |
---|
| 183 | ! run standard (ispinup=0) et pour les initialisations |
---|
| 184 | !--------------------------------------------------------------------------------- |
---|
| 185 | |
---|
| 186 | if (itracebug.eq.1) call tracebug(' Dans spinup_3_bed') |
---|
| 187 | |
---|
| 188 | if ((nbed.eq.1).and.nt.ne.1.and.isynchro.eq.1) then |
---|
| 189 | call bedrock() ! bedrock adjustment |
---|
| 190 | endif |
---|
| 191 | end if spinup_3_bed |
---|
| 192 | |
---|
| 193 | spinup_4_vitdyn: if (ispinup.ge.1.or.nt.lt.5.) then ATTENTION A MODIFIER ipsinup=3 shelf_vitbil |
---|
| 194 | |
---|
| 195 | !--------------------------------------------------------------------------------- |
---|
| 196 | ! spinup_4_vitdyn |
---|
| 197 | !---------------- |
---|
| 198 | ! Ne passe dans ce bloc que quand on veut calculer les vitesses dynamiques : |
---|
| 199 | ! - initialisations (nt.lt.5) |
---|
| 200 | ! - run standard (ispinup=0) (y compris iter_beta.ne.0) |
---|
| 201 | ! - spinup temperature avec vitesses dynamiques (ispinup = 1) |
---|
| 202 | !--------------------------------------------------------------------------------- |
---|
| 203 | |
---|
| 204 | ! (*********** DIFFUSIVITIES - VITESSES SHELF **********************) |
---|
| 205 | |
---|
| 206 | |
---|
| 207 | call Neffect() |
---|
| 208 | |
---|
| 209 | call diffusiv() |
---|
| 210 | |
---|
| 211 | |
---|
| 212 | iterbeta: if (iter_beta.ne.0) then ! pour les iterations sur beta |
---|
| 213 | |
---|
| 214 | if (.not.shelf_vitbil) then ! bloc si les vitesses shelves viennent de diagnoshelf |
---|
| 215 | |
---|
| 216 | ! Vcol declare dans le module spinup_mod |
---|
| 217 | ! ou dans le dragging si no_spinup |
---|
| 218 | |
---|
| 219 | where (flotmx(:,:)) |
---|
| 220 | uxbar(:,:)=uxflgz(:,:) |
---|
| 221 | elsewhere |
---|
| 222 | uxbar(:,:)=Vcol_x(:,:) |
---|
| 223 | end where |
---|
| 224 | |
---|
| 225 | where (flotmy(:,:)) |
---|
| 226 | uybar(:,:)=uyflgz(:,:) |
---|
| 227 | elsewhere |
---|
| 228 | uybar(:,:)=Vcol_y(:,:) |
---|
| 229 | end where |
---|
| 230 | else |
---|
| 231 | uxbar(:,:)=Vcol_x(:,:) |
---|
| 232 | uybar(:,:)=Vcol_y(:,:) |
---|
| 233 | end if |
---|
| 234 | end if iterbeta |
---|
| 235 | |
---|
| 236 | |
---|
| 237 | call strain_rate() |
---|
| 238 | |
---|
| 239 | if (iter_beta.ne.0) then |
---|
| 240 | uxbar(:,:)=uxflgz(:,:)+uxdef(:,:) |
---|
| 241 | uybar(:,:)=uyflgz(:,:)+uydef(:,:) |
---|
| 242 | end if |
---|
| 243 | |
---|
| 244 | |
---|
| 245 | |
---|
| 246 | !----------------------------------------------------------------------- |
---|
| 247 | ! debut bloc appels tous les dtt |
---|
| 248 | !------------- ------------------------------------------------------- |
---|
| 249 | |
---|
| 250 | if ((shelfy).and.(marine).and.(isynchro.eq.1)) then |
---|
| 251 | |
---|
| 252 | call lake_to_slv |
---|
| 253 | !write(6,*) 'premier diagno dans step' |
---|
| 254 | !call flottab |
---|
| 255 | |
---|
| 256 | ! les 3 lignes ci-dessous pour avoir les champs avant diagno (en cas de pb colineaire) |
---|
| 257 | !call init_sortie_ncdf |
---|
| 258 | !call sortie_ncdf_cat |
---|
| 259 | |
---|
| 260 | call diagnoshelf |
---|
| 261 | |
---|
| 262 | !call sortie_ncdf_cat |
---|
| 263 | ! STOP |
---|
| 264 | |
---|
| 265 | |
---|
| 266 | ! faire 10 tours pour equilibrer |
---|
| 267 | if ((icompteur.eq.0).and.(nt.lt.10)) then |
---|
| 268 | do m=1,1 |
---|
| 269 | ! write(6,*) 'dans boucle m' |
---|
| 270 | call diagnoshelf() |
---|
| 271 | end do |
---|
| 272 | end if |
---|
| 273 | |
---|
| 274 | endif ! |
---|
| 275 | |
---|
| 276 | call mix_SIA_L1 |
---|
| 277 | |
---|
| 278 | debug_3D(:,:,78)= abs(uxbar(:,:))-abs(Vcol_x(:,:)) |
---|
| 279 | debug_3D(:,:,79)= abs(uybar(:,:))-abs(Vcol_y(:,:)) |
---|
| 280 | |
---|
| 281 | debug_3D(:,:,80)= abs(uxdef(:,:))-abs(Vcol_x(:,:)) |
---|
| 282 | debug_3D(:,:,81)= abs(uydef(:,:))-abs(Vcol_y(:,:)) |
---|
| 283 | |
---|
| 284 | |
---|
| 285 | end if spinup_3_bed |
---|
| 286 | |
---|
| 287 | end subroutine step_grisli |
---|
| 288 | |
---|
| 289 | !-------------------------------------------------------------------------- |
---|
| 290 | ! step_grisli1 est a l'interieur de la boucle temps. |
---|
| 291 | ! il commence par icethick, qui détermine dt, isynchro |
---|
| 292 | |
---|
| 293 | ! ensuite contient les appels aux diverses sorties |
---|
| 294 | ! puis un appel a step_grisli |
---|
| 295 | |
---|
| 296 | |
---|
| 297 | subroutine step_grisli1 |
---|
| 298 | |
---|
| 299 | use module3d_phy |
---|
| 300 | |
---|
| 301 | use geom_type |
---|
| 302 | use temperature_type |
---|
| 303 | use ice_flow_type |
---|
| 304 | use mask_flgz_type |
---|
| 305 | use deformation_type |
---|
| 306 | use autre_pr_temp_type |
---|
| 307 | |
---|
| 308 | use module_choix ! module de choix du type de run |
---|
| 309 | ! module_choix donne acces a tous les modules |
---|
| 310 | ! de declaration des packages |
---|
| 311 | use sorties_ncdf_grisli |
---|
| 312 | use flottab_mod |
---|
| 313 | use interface_icetempmod |
---|
| 314 | use diagno_mod |
---|
| 315 | use track_debug |
---|
| 316 | |
---|
| 317 | if (itracebug.eq.1) write(num_tracebug,*) ' Entree dans step_grisli1', & |
---|
| 318 | ' nt=',nt,'iter_beta=',iter_beta |
---|
| 319 | |
---|
| 320 | nospinup_4: if ((ispinup.eq.0.or.nt.lt.5).and.(iter_beta.eq.0)) then |
---|
| 321 | if (itracebug.eq.1) call tracebug(' Dans nospinup_4 : appelle icethick') |
---|
| 322 | |
---|
| 323 | call icethick3() |
---|
| 324 | |
---|
| 325 | call calving |
---|
| 326 | |
---|
| 327 | if (itracebug.eq.1) call tracebug('apres calving') |
---|
| 328 | |
---|
| 329 | else if ((ispinup.eq.1).or.(ispinup.eq.3).or.(iter_beta.gt.0)) then |
---|
| 330 | dt=dtmax |
---|
| 331 | hdot(:,:)=0. |
---|
| 332 | |
---|
| 333 | ! if (nt.eq.6) time=0. |
---|
| 334 | ! time=nint(dtmax*nint(time/dtmax))+dtmax semble poser des problemes |
---|
| 335 | time=nt*dtmax |
---|
| 336 | |
---|
| 337 | if (itracebug.eq.1) write(num_tracebug,*) ' Dans nospinup_4 no icethick ',ispinup, & |
---|
| 338 | ' time=',time |
---|
| 339 | else if (ispinup.eq.2) then |
---|
| 340 | call force_balance_vel |
---|
| 341 | call icethick3() |
---|
| 342 | call calving |
---|
| 343 | |
---|
| 344 | end if nospinup_4 |
---|
| 345 | |
---|
| 346 | CALL write_trace_940('ICE SHEET SURFACE') |
---|
| 347 | |
---|
| 348 | where(.not.flot(:,:)) |
---|
| 349 | S(:,:)=Bsoc(:,:)+H(:,:) |
---|
| 350 | B(:,:)=Bsoc(:,:) |
---|
| 351 | end where |
---|
| 352 | |
---|
| 353 | ! calcul de Hmx et Hmy -> shift=-1, dim=1 -> H(i-1,j) |
---|
| 354 | |
---|
| 355 | hmx(:,:)=0.5*(H(:,:)+eoshift(H(:,:),shift=-1,boundary=0.,dim=1)) |
---|
| 356 | hmy(:,:)=0.5*(H(:,:)+eoshift(H(:,:),shift=-1,boundary=0.,dim=2)) |
---|
| 357 | hmx(1,:)=0. |
---|
| 358 | hmy(:,1)=0. |
---|
| 359 | |
---|
| 360 | ! mise a 0 des vitesses quand epaisseur nulle |
---|
| 361 | where (hmx(:,:).le.1.) |
---|
| 362 | uxbar(:,:)=0. |
---|
| 363 | endwhere |
---|
| 364 | where (hmy(:,:).le.1.) |
---|
| 365 | uybar(:,:)=0. |
---|
| 366 | endwhere |
---|
| 367 | |
---|
| 368 | ! ________________________________ fin bloc appel tous les tend |
---|
| 369 | |
---|
| 370 | ! (****** ICE SHEET STATISTICS and SHORT OUTPUT *******) |
---|
| 371 | ! if (mod(nt,nint(NDISP/DT)).eq.0).or.(nt.eq.1).or.(nt.eq.2).or. |
---|
| 372 | ! & (nt.eq.3).or.(nt.eq.4).or.(nt.eq.5).or.(nt.eq.10).or. |
---|
| 373 | |
---|
| 374 | ! if ((mod(abs(TIME),NDISP*1.).lt.dtmin).or.(NT.le.5).or.(NT.eq.10) & |
---|
| 375 | ! .or.(nt.eq.15).or.(nt.eq.20).or.(nt.eq.30).or. & |
---|
| 376 | ! (abs(TIME-TEND).lt.dtmin)) then |
---|
| 377 | |
---|
| 378 | ! call sortie_ncdf_cat |
---|
| 379 | |
---|
| 380 | |
---|
| 381 | if ((mod(abs(TIME),NDISP*1.).lt.dtmin).or.(isynchro.eq.1)) then |
---|
| 382 | |
---|
| 383 | call shortoutput() |
---|
| 384 | |
---|
| 385 | if ((geoplace.eq.'anteis1') & |
---|
| 386 | .or.(geoplace.eq.'ant20km')) then |
---|
| 387 | |
---|
| 388 | call ts_out() |
---|
| 389 | endif |
---|
| 390 | ! call sealevel_out |
---|
| 391 | |
---|
| 392 | endif |
---|
| 393 | |
---|
| 394 | |
---|
| 395 | ! Sorties horizontales. |
---|
| 396 | !-------------------------------------------------------------------------------- |
---|
| 397 | ! Premiere partie des sorties horizontales |
---|
| 398 | |
---|
| 399 | ! if (mod(abs(dble(TIME)),dble(DTSORTIE)).lt.dble(dtmin)) then |
---|
| 400 | ! |
---|
| 401 | |
---|
| 402 | if ((mod(abs(TIME),dtt).lt.dtmin).or.(isynchro.eq.1)) then |
---|
| 403 | isynchro=1 |
---|
| 404 | ! if (isynchr o.eq.1) then |
---|
| 405 | call testsort_time(real(time)) |
---|
| 406 | if (iglob_hz.eq.1) then |
---|
| 407 | call sortie_hz_multi ! pour les variables déclarées dans 3D |
---|
| 408 | ! call hz_output(real(time)) |
---|
| 409 | endif |
---|
| 410 | |
---|
| 411 | else |
---|
| 412 | iglob_hz=0 |
---|
| 413 | |
---|
| 414 | endif |
---|
| 415 | |
---|
| 416 | |
---|
| 417 | ! *** Sortie des profils tous les dtprofile *** |
---|
| 418 | if ((NT.eq.1) & ! .or.(NT.eq.2) |
---|
| 419 | .or.(mod(abs(dble(TIME)),dble(DTPROFILE)).lt.dble(dtmin)) & |
---|
| 420 | .or.(abs(TIME-TEND).lt.dtmin)) then |
---|
| 421 | |
---|
| 422 | ! call sortieprofile() |
---|
| 423 | |
---|
| 424 | endif |
---|
| 425 | |
---|
| 426 | ! *** SORTIE COMPTEUR TOUS LES DTCPT ANS *** |
---|
| 427 | |
---|
| 428 | if ((mod(abs(dble(TIME)),dble(DTCPT)).lt.dble(dtmin)).OR. & |
---|
| 429 | (ABS(TIME+297000).LT.dtmin).OR. & |
---|
| 430 | (ABS(TIME+126000).LT.dtmin).OR. & |
---|
| 431 | (ABS(TIME+21000).LT.dtmin).OR. & |
---|
| 432 | (ABS(TIME+16000).LT.dtmin)) THEN |
---|
| 433 | |
---|
| 434 | call sortiecptr() |
---|
| 435 | |
---|
| 436 | endif |
---|
| 437 | |
---|
| 438 | |
---|
| 439 | ! ************* sorties netcdf **************** |
---|
| 440 | |
---|
| 441 | call testsort_time_ncdf(real(time)) |
---|
| 442 | if (iglob_ncdf .EQ. 1) call sortie_ncdf_cat |
---|
| 443 | |
---|
| 444 | ! call sortie_ncdf_cat |
---|
| 445 | |
---|
| 446 | call track_change_T |
---|
| 447 | |
---|
| 448 | call step_grisli() |
---|
| 449 | ! |
---|
| 450 | ! |
---|
| 451 | ! |
---|
| 452 | end subroutine step_grisli1 |
---|