| 1 | subroutine spreadglacier_simple(ngrid,nq,pqsurf,timstep) |
|---|
| 2 | |
|---|
| 3 | use geometry_mod, only: latitude, longitude, cell_area |
|---|
| 4 | use radinc_h, only : naerkind |
|---|
| 5 | use surfdat_h, only: phisfi |
|---|
| 6 | use tracer_h, only: igcm_n2 |
|---|
| 7 | |
|---|
| 8 | implicit none |
|---|
| 9 | |
|---|
| 10 | !================================================================== |
|---|
| 11 | ! Purpose |
|---|
| 12 | ! ------- |
|---|
| 13 | ! Spreading the glacier in a circular crater using a speed of flow |
|---|
| 14 | ! |
|---|
| 15 | ! Inputs |
|---|
| 16 | ! ------ |
|---|
| 17 | ! ngrid Number of vertical columns |
|---|
| 18 | ! pqsurf(ngrid,nq) N2 ice tracer on surface kg/m2 |
|---|
| 19 | ! |
|---|
| 20 | ! Outputs |
|---|
| 21 | ! ------- |
|---|
| 22 | ! pqsurf(ngrid,nq) new value for N2 ice tracer on surface kg/m2 |
|---|
| 23 | ! |
|---|
| 24 | ! Authors |
|---|
| 25 | ! ------- |
|---|
| 26 | ! Tanguy Bertrand (2015) |
|---|
| 27 | ! |
|---|
| 28 | !================================================================== |
|---|
| 29 | |
|---|
| 30 | #include "dimensions.h" |
|---|
| 31 | |
|---|
| 32 | !----------------------------------------------------------------------- |
|---|
| 33 | ! Arguments |
|---|
| 34 | |
|---|
| 35 | INTEGER ngrid, nq, ig, ig2 |
|---|
| 36 | REAL :: pqsurf(ngrid,nq) |
|---|
| 37 | REAL timstep |
|---|
| 38 | !----------------------------------------------------------------------- |
|---|
| 39 | ! Local variables |
|---|
| 40 | |
|---|
| 41 | REAL :: zdqsurf(ngrid) |
|---|
| 42 | REAL phiref |
|---|
| 43 | REAL dist |
|---|
| 44 | REAL locdist |
|---|
| 45 | REAL tau |
|---|
| 46 | REAL speed |
|---|
| 47 | REAL dist_pluto2 ! function |
|---|
| 48 | LOGICAL firstcall |
|---|
| 49 | DATA firstcall/.true./ |
|---|
| 50 | INTEGER,SAVE,ALLOCATABLE :: indglac(:) |
|---|
| 51 | INTEGER, SAVE :: nglac |
|---|
| 52 | |
|---|
| 53 | !---------------- INPUT ------------------------------------------------ |
|---|
| 54 | ! Defined for Tombaugh crater : |
|---|
| 55 | |
|---|
| 56 | dist=300. ! reference radius km |
|---|
| 57 | phiref=-1000. ! geopotential reference |
|---|
| 58 | ! tau=1.e10 |
|---|
| 59 | speed=6.4e-7 ! glacier speed km/s |
|---|
| 60 | ! tau=1. ! direct redistribution glacier |
|---|
| 61 | |
|---|
| 62 | !--------------- FIRSTCALL : define crater index ----------------------- |
|---|
| 63 | IF (firstcall) then |
|---|
| 64 | ALLOCATE(indglac(ngrid)) |
|---|
| 65 | nglac=0 |
|---|
| 66 | indglac(:)=0 |
|---|
| 67 | DO ig=1,ngrid |
|---|
| 68 | IF (phisfi(ig).lt.phiref) THEN |
|---|
| 69 | nglac=nglac+1 |
|---|
| 70 | indglac(nglac)=ig |
|---|
| 71 | ENDIF |
|---|
| 72 | ENDDO |
|---|
| 73 | firstcall=.false. |
|---|
| 74 | |
|---|
| 75 | ENDIF |
|---|
| 76 | !--------------- redistribution ------------------------------------- |
|---|
| 77 | DO ig=1,nglac |
|---|
| 78 | ! redistribution if large amount of ice only : detection of reservoirs |
|---|
| 79 | IF (pqsurf(indglac(ig),igcm_n2).gt.50.) THEN |
|---|
| 80 | ! looking where to redistribute : if less ice, lower altitude |
|---|
| 81 | ! and close to the reservoir |
|---|
| 82 | DO ig2=1,nglac |
|---|
| 83 | IF ((pqsurf(indglac(ig2),igcm_n2).lt.pqsurf(indglac(ig),igcm_n2)*0.8)) THEN |
|---|
| 84 | !write(*,*) 'loop=',latitude(indglac(ig2))*180/3.14,longitude(indglac(ig2))*180/3.14 |
|---|
| 85 | !.and.(phisfi(indglac(ig2)).le.phisfi(indglac(ig)))) THEN |
|---|
| 86 | locdist=dist_pluto2(latitude(indglac(ig2)),longitude(indglac(ig2)),latitude(indglac(ig)),longitude(indglac(ig))) |
|---|
| 87 | !write(*,*) 'dist=',locdist |
|---|
| 88 | IF (locdist.lt.dist) THEN |
|---|
| 89 | !write(*,*) 'calcul=',latitude(indglac(ig2))*180/3.14,longitude(indglac(ig2))*180/3.14,latitude(indglac(ig))*180/3.14,longitude(indglac(ig))*180/3.14 |
|---|
| 90 | tau=locdist/speed ! characteristic time for redistribution |
|---|
| 91 | zdqsurf(indglac(ig2))=pqsurf(indglac(ig),igcm_n2)*(1.-exp(-timstep/tau))/timstep*cell_area(indglac(ig))/cell_area(indglac(ig2)) |
|---|
| 92 | zdqsurf(indglac(ig))=pqsurf(indglac(ig),igcm_n2)*(exp(-timstep/tau)-1.)/timstep |
|---|
| 93 | pqsurf(indglac(ig2),igcm_n2)=pqsurf(indglac(ig2),igcm_n2)+zdqsurf(indglac(ig2))*timstep |
|---|
| 94 | pqsurf(indglac(ig),igcm_n2)=pqsurf(indglac(ig),igcm_n2)+zdqsurf(indglac(ig))*timstep |
|---|
| 95 | ENDIF |
|---|
| 96 | ENDIF |
|---|
| 97 | ENDDO |
|---|
| 98 | ENDIF |
|---|
| 99 | ENDDO |
|---|
| 100 | |
|---|
| 101 | end subroutine spreadglacier_simple |
|---|
| 102 | |
|---|
| 103 | real function dist_pluto2(lat1,lon1,lat2,lon2) |
|---|
| 104 | implicit none |
|---|
| 105 | !#include "comcstfi.h" |
|---|
| 106 | real, intent(in) :: lat1,lon1,lat2,lon2 |
|---|
| 107 | real dlon,dlat |
|---|
| 108 | real a,c,d |
|---|
| 109 | real rad |
|---|
| 110 | rad=1187. ! planetary radius (km) 1187 km sur Pluton |
|---|
| 111 | |
|---|
| 112 | dlon = lon2 - lon1 |
|---|
| 113 | dlat = lat2 - lat1 |
|---|
| 114 | a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2 |
|---|
| 115 | c = 2 * atan2( sqrt(a), sqrt(1-a) ) |
|---|
| 116 | dist_pluto2 = rad * c |
|---|
| 117 | return |
|---|
| 118 | end function dist_pluto2 |
|---|
| 119 | |
|---|
| 120 | |
|---|
