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