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